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The vartous tasks described in thus report dealt with numerical and 
experimental analysis of a thin liquid film on a rotating and statronary 
disk related to the development of an absorber unit for a high capacrty 
spacecraft absorptron coolrng system. The idea that was focused upon tn 

this report was the creation of an artificial gravity by the use of a 

centrifugal field. The basic phenomena related to the fluid flow and 

heat transfer on rotating systems that have been reported during this 

effort can be applied to other areas of space systems as well. 
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Section I 


ORE- DDffiESIORAL COIPUTATIOR OP THIN LIQUID FILI FLOWS 


1.1 S Till A RY 


The flow of a thin liquid film with a free surface along a horizontal 
plate which emanates from a pressurized vessel is examined numerically. In 
one g, a hydraulic jump was predicted in both plane and radial flow, which 
could be forced away from the inlet by increasing the inlet Froude number 
or Reynolds number. In zero g, the hydraulic jump was not predicted. The 
effect of solid- body rotation for radial flow in one g was to "wash out 
the hydraulic jump and to decrease the film height on the disk. The liquid 
film height under one g and zero g was equal under solid- body rotation 
because the effect of centrifugal force was much greater than that of the 
gravitational force. The heat transfer to a film on a rotating disk was 
predicted to be greater than that of a stationary disk because the liquid 
film is extremely thin and is moving with a very high velocity. 
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1.2 INTRODUCTION 


The study of thin liquid films has been performed by many researchers 
in the past due to the high heat transfer rates that can be achieved. Even 
though the heat transfer to thin films falling down a vertical plate is 
quite high, much better heat transfer rates can be realized when the film 
is generated by the impingement of a liquid jet onto a surface that is 
perpendicular to the jet. This is due to the fact that in most cases the 
film velocities are greater than those of a falling film. Since it is 
difficult to examine the entire flow field from the impinging jet to the 
perpendicular thin film, it is proposed to study only the thin film. This 
is accomplished by the following mechanism: the liquid is pressurized in a 
container and emanates from a slot which generates the thin film on a 
horizontal plate. Vith this situation, the height and mean velocity of the 
film will be known boundary conditions, which will aid in determining the 
characteristics of the film downstream. After the structure of the liquid 
film has been examined, the amount of heat that can be transferred to it 
can be found. Furthermore, in this study it is also desired to understand 
the effects of centrifugal force and zero gravity on the flow field. If 
the film is generated on a spinning disk, it is expected that the film will 
accelerate and become thinner, since the centrifugal force acts in the main 
direction of the flow. This acceleration of the liquid film will increase 
the amount of heat that can be transferred from the spinning disk to the 
film. The effect of zero g on the flow is also examined to determine the 
nature of the liquid film height when the gravitational field across the 
film is not present. The combined effects of centrifugal acceleration and 
a zero g environment on the fluid characteristics and heat transfer to a 


2 



thin film are not known at present. It is felt, however, that this 
fundamental research will become important in the future because of 
programs like the Space Station where improved heat exchangers will be 
needed. This research will also be instrumental in the development of an 
absorber unit for a spacecraft vapor- absorption heat pump system, which was 
the motivation for the present work. 

Many investigations have been carried out concerning thin liquid films 
created by impinging jets and/or spinning disks because of the numerous 
applications where the improved .heat or mass transfer is desirable. 
Sparrow and Gregg (1959) analytically determined the condensate thickness 
on a rotating disk in a large vessel of quiescent saturated vapor. Vatson 
(1964) analysed a free- falling jet which impinges on a horizontal plane 
using a similarity solution. Chaudhury (1964) used a similarity solution 
to analyse the heat transfer in the thin liquid layer on a horizontal plate 
on which a liquid jet impinges. Aroesty et al. (1967) studied the use of 
thin films under centrifugal force as an aid to blood oxygenation. The 
flow of a liquid film on a rotating disk was examined by Matsumoto et al. 
(1973) in connection with the atomization of a liquid for spray drying and 
the promotion of chemical reactions or absorption between gases and 
liquids. Rauscher et al. (1973) analysed the laminar flow of a thin film 
on a rotating disk by employing the asymptotic expansion technique. 
Eliseev (1983) examined the spatial stability of liquid films on a rotating 
disk by the asymptotic expansion technique. The study of thin films in 
regard to spin coating deposition of thin solid films was carried out by 
Jenekhe (1984). Needham and Merkin (1987) theoretically studied thin 
axisymmetric liquid films on a horizontally rotating disk to determine the 
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criteria for stability of the film. Katto and Yokoya (1988) analyzed the 

ft 

existing experimental data for the critical heat flux of a disk heater 
cooled by an impinging liquid jet and gave an equation correlating the 
data. 


The motivation of the present work is the study of thin liquid films 
which will be used in space- based centrifugal heat exchangers and vapor 
absorption heat pumps. Since a falling film cannot be achieved in a zero g 
environment, it is natural to consider the possibility of using the 
artificial gravity created by the. centrifugal force on a rotating disk to 
generate a thin film. 

In the present study, the governing equations and boundary conditions 
are presented for two situations of a thin liquid layer emanating from a 
pressurized vessel and traveling along a horizontal plate with a constant 
initial height and uniform initial velocity as shown in Fig. 1.1. The 
first case is when the liquid flows along a channel with a constant width, 
which is shown in Fig. 1.2a. This situation is the same as open channel 
flow, but since the liquid height is very thin the effect of viscosity must 
be accounted for. The second case, which is shown in Fig. 1.2b, is when 
the liquid originates between two parallel disks, and then spreads out 
radially over the bottom disk in a free surface thin film because the 
diameter of the upper disk is much smaller than that of the lower disk. 
This situation is similar to a jet of liquid impinging onto the center of a 
horizontal disk, except that the height and the mean velocity of the liquid 
film are known at a specific radial location. In the case of radial flow, 
the rotational acceleration will be modeled as solid- body rotation, where 
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the centrifugal forces act only in the radial ^direction. 

It is desired to numerically solve for the liquid height at any 
distance down the length of the plate for different Froude numbers and 
Reynolds numbers specified at the inlet. The heat transfer from the plate 
will also be studied for the case of simple heating with no evaporation at 
the free surface. Since the inlet Froude number may be greater than unity, 
it is possible that a hydraulic jump will occur at some point in the 
computational domain. A hydraulic jump is when the flow suddenly changes 
from supercritical (Fr > 1) to. subcritical (Fr < 1) flow, which is 
accompanied by a sudden increase in the liquid height. This is analogous 
to the shock wave in gas dynamics when the flow changes from supersonic 
(M > 1) to subsonic (M < 1) flow in a very short distance. The similarity 
between the hydraulic jump and the shock wave in gas dynamics suggests 
using the familiar approach of modeling the flow as a transient phenomenon 
and allowing the solution to march in time to achieve the desired 
steady- state results. To the authors’ knowledge, this approach has not 
been applied to this type of problem before. It should be emphasized that 
conventional numerical techniques that are applied to 2- or 3- dimensional 
problems cannot be used directly to solve this particular problem because 
the location of the free surface is unknown prior to the start of the 
calculations . 

1.3 FLDID FLQV ANALYSIS 

The conservation of mass, momentum, and energy in the general 
coordinate system can be expressed in the following compact form for a thin 
incompressible liquid film with constant properties. 


7 


7 • V = 0 


(1.3.1) 


4 = - Vp ♦ ,V 2 V ♦ 


P § 


(1.3.2) 


DT 

Dt 


= aV 2 T 


(1.3.3) 


In the general coordinate system given in Fig. 1.1, the direction that 
is parallel to the plate in the main direction of the flow will be denoted 
by r and the direction that is normal to the plate will be denoted by z. 

The velocity vector V has three components; w, v, and u, where w is in the 
r- direction and v is in the y- direction. The component of velocity in the 
x-direction, as shown in Fig. 1.2a, is set equal to zero for plane flow. 
The component of velocity in the 0- direction, as shown in Fig. 1.2b, is 
equal to u for radial flow. 


Several assumptions are made in order to reduce the complexity of the 
governing equations. Any changes in the variables with respect to the 9 
or x-directions are set to zero. In the gravitational field, g, only the 
component across the film thickness, g^, is considered. The boundary layer 
assumptions are imposed on the governing equations, which neglect the 
diffusion terms in the main direction of the flow based on an 
order- of- magnitude analysis. At the free surface, the effects of 
interfacial shear stress and the surface tension are assumed to be zero. 
It should be noted that the effect of surface tension may be significant 
near the inlet and outlet of the domain and in the hydraulic jump region. 
The pressure in the liquid film is defined as the difference between the 
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actual pressure and the atmospheric pressure, so at the free surface the 
pressure is equal to zero. For the case of radial flow, it is assumed that 
the velocity component in the ^-direction is u = wr, i.e., solid-body 
rotation. By an order- of- magnitude analysis, this assumption is valid for 
small values of the Rossby number (Ro < 1). 

The boundary and initial conditions for the governing equations in 
fixed coordinates are those on the flat plate and on the free surface of 
the liquid. 


y = 0: 

w = 0, 

u = wr, 

v = 0 


86 

86 _ 

II 

O 

y = 6: 

v = 3t + w 3r’ P = 

r = r lS 

w w in 5 

v = 0, 

f ■ *i„ 

t = 0: 

W = w 0 , V = 0 



The boundary condition involving v at the free surface is the 
kinematic condition for time- dependent problems. The boundary condition 
involving dvjdy at y = 6 corresponds to neglecting the interfacial shear 
stress at the free surface. 


The continuity equation and the conservation of momentum equations in 
the r- and y- directions can be integrated with respect to the y- direction 
across the thin liquid layer by using Leibniz’s rule and the kinematic 
condition at the free surface. Due to the boundary- layer assumptions, the 
conservation of momentum equation in the y- direction results in a balance 
between the pressure gradient and the gravitational force. This equation 
is then substituted into the r- direction momentum equation. It is assumed 


9 



that the velocity in the main direction of the flow- is constant across the 
thin film, so that 

(6 2 2 
w J dy = V6 
J o 


It should be noted that this is not a major assumption. For example, 

2 

for a parabolic velocity profile, the left-hand side is equal to 6V 6/5. 
The resulting one- dimensional time- dependent equations in the general 
coordinate system are given in the following form: 

T I? * M - 0 (1 ' 3 ' 4) 

r 


d 

3t 


W) 


d 

<?r 


V 2 6 


1 *2 
^6 


\6 


r 2 

ra; 



r 


1 <9v 


y = 0 


(1.3.5) 


To properly assign the boundary conditions of the flow, the 
characteristic behavior of the flow will be examined. The governing 
equations can be cast in the characteristic form by the following analysis. 


The governing equations (1.3.4 and 1.3.5) can be rearranged by 
subtracting V continuity equation from the momentum equation. 


1 36 V 36 3V _ IV 

7H + 7^r + ^r = ' T 


3M 

M + 



(g<5) 


1 36 
7 Jr 


= TV 


r 

v 
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The characteristic form of the equations are obtained by linearizing 
the continuity and momentum equations: 


(It + V inl?]l“ * l? (V) = S 1 


( 


d , v d 

3t “in?? 


u .2 3 / S'. 
» — W (l-> 


in 


= S r 


where 


a. = g 6- 
m 6 in 

q _ AV 

S 1 " ’ r 

s 2 = rw 2 [l - 


Spru* 


These two equations can be represented by the following characteristic 
equation. 


<?R. dl. 

W + c i W ~ S i 


' The Reimann invariants and the physical wave speeds are given in the 
following table. 


v * a in W* in ) 
V - » ta (*/*i.) 


c i 

V. ♦ a. 
m m 

V. - a. 
in m 


S i 

a in S l + S 2 
a in S l + S 2 
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An examination of the characteristics reveals that the first invariant 
will always propagate downstream. The direction of propagation for the 
second invariant depends upon whether the flow is supercritical or 
subcritical. For supercritical flow (i.e., Fr > 1), the second invariant 
propagates downstream, so the boundary conditions on V and 6 will both be 
assigned at the inlet of the flow field. For subcritical flow (i.e. , 
Fr < 1), the second invariant propagates upstream, so the boundary 
conditions must be assigned at the inlet and the outlet of the flow field. 


The following assumption is made to simplify the governing equations 
so that the film height can be eliminated from the equations. 



0 


This assumption results in the following quasi - steady governing 
equations : 


l X MS = 


(2r) 


^ j — constant 


rdV 

5 tt 


d 

3F 


\v 2 6 + 


a*' 4 *! 


= xs 


TU 



( 1 . 3 . 6 ) 


( 1 . 3 . 7 ) 


By using this assumption, it is understood that the unsteady solutions 
with respect to time are not precisely accurate, but the accuracy of the 
steady- state solution is not affected. 
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The following analysis for the skin- friction coefficient is given 
since the shear stress at the wall is not known. When the boundary layer 
on the plate is much smaller than the film thickness, i.e., developing 
flow, it is assumed that the boundary layer is the same as when there is a 
uniform far- field boundary. Therefore, when 6 > <5g, where <5g is the 
Blasius boundary layer thickness, the Blasius skin friction coefficient 
will be employed to approximate the shear stress term. Vhen the boundary 
layer reaches the free surface (i.e., 6 = £g) , a parabolic velocity profile 
is used to approximate the shear stress term. The Blasius boundary layer 
thickness is given by 



5.0 



The Blasius skin friction coefficient is as follows: 


0.664 



The Blasius skin friction coefficient can also be presented in terms 
of the Blasius boundary layer thickness. 

/ n 1 \ 3 . 32 

Cr(Blasius) = 

f^B' 

v J 

For 6 = it is assumed that the velocity profile is parabolic. 
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The coefficient of skin friction is 


^(parabolic) 



dv 

% 
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The two skin friction coefficients match when 


c f (Blasius) = c f (parabolic) 


or 



1.81 


The value of the skin friction will be determined in the following manner. 


5 -°( r - r J 

V ( r • r in) 


0 < 


6 

< 1.81 



5 


> 


1.81 


0.664 



The governing equations are nondimensionalized in the following manner 


to simplify the analysis and to generalize the results. 

V. ' 


V_ 

V. 


= V 


in 


5 r + 

TT = 1 

in 




= < 


^ m 

= r 

in 


V 9 

— i— = Fr 2 

PT7 in 


vj. 

1 m 


= Re. 
in 


in 


V. 

-12- = Ro- 
ur . m 

in 


On 

_5_ - = dimensionless boundary layer thickness 

o ■ u 

in 


The Reynolds number of the flow was chosen to nondimensionalize the 
governing equations instead of the Ekman number because the Reynolds number 
is needed when the flow is stationary. 


After eliminating the dimensionless film height from the momentum 
equation with the continuity equation, the dimensionless governing equation 
is as follows: 


v + v^g; = H 

T S t 


( 1 . 3 . 9 ) 


where 


G 



.ypJ 


/n; 


( 1 . 3 . 10 ) 


and 




U 

{ 


for 


0 < — < 1.81 


( 1 . 3 . 11 ) 
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or 


H = A 




*!] . 0-332 V 5/2 |j_ f for S - > 1.81 (1.3.12) 

< J ^in $ 


where the dimensionless boundary layer thickness is given by 


= 5-0 


« - <L> 


"VTEeT 


in 


(1.3.13) 


The dimensionless initial and boundary conditions are given as follows: 


▼(o,b = i 


V( T ^ in ) = 1 


These conditions correspond to a constant velocity across the region at the 
initial time and a constant inlet velocity at all times. Since the 
steady- state results are the only concern due to the constraints made on 
the mathematical modeling (36/ dt ~ 0), the specification of the initial 
condition is required only for the numerical iteration scheme. 


An artificial viscosity term (/JV^) is included in the governing 
equation to dampen numerical oscillations in the solution (Anderson et al., 
1984). The oscillations are commonly referred to as the Gibb’s phenomenon, 
which occurs near a double- valued point in the solution, such as a 
hydraulic jump or a shock wave. This term reduces the numerical 
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oscillations in the immediate vicinity of the hydraulic jump, but does not 
affect the numerical solution elsewhere in the domain. 

1.4 HEAT TRANSFER ANALYSIS 


It is desired to find the amount of heat transferred from the plate or 
disk to the liquid when the plate is heated. This analysis assumes that 
the velocity of the thin film is approximated by the similarity profile 
presented in equation (1.3.8). The temperature profile across the thin 
film is assumed to be a quadratic function, with the coefficients 
determined by the boundary conditions imposed at the surface of the plate 
and the free surface. The boundary conditions are a constant heat flux at 
the plate and the free surface of the liquid film is adiabatic, which 
corresponds to simple heating with no evaporation at the free surface. The 
boundary conditions are given as: 


dfl 

3y 


y = 6 


= 0 


(1.4.1) 


-K 


dT 

Jy 


y = 0 


(1.4.2) 


The temperature distribution across the thin film is 

i-Vt'fsb * 2 (1 - 4 - 3) 

The mixed- mean temperature of the film is given by the following equation 



(1.4.4) 


T - ^ 
l b ' TTC 





The heat transfer coefficient is defined by 



5 K 
2 7 


(1.4.5) 


Therefore, the heat transfer coefficient can be found in terms of the 
liquid film height. 


The modified Nusselt number is defined in the following manner: 




1/3 


2.5 

T" 


2 , 


gj 


1/3 


(1.4.6) 


This definition has been used previously in the literature concerning thin 
falling films. For a zero g situation, the Nusselt number is defined, in 
terms of the liquid film thickness: 


Nu*=^=2.5 (1-4.7) 

For practical applications, the mixed- mean temperature of the fluid should 
be known. Therefore, the following analysis is given to find the 
mixed- mean temperature for the case of simple heating with no evaporation 
at the free surface for plane and radial flow. This analysis follows the 
same pattern as was accomplished for the momentum equation of the problem. 
The conservation of energy will be integrated across the thin film in order 
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to find the mixed- mean temperature. 


The conservation of energy equation is 


P c r 


m 


31 9T 

K d 

A0T] 

+ W 3F + 

= 

r 3r 


+ K 


d 2 1 

T2 

dy 


( 1 . 4 . 8 ) 


The assumption is made that the diffusive term in the r- direction is zero. 


5T 


+ 



+ v 


5T _ 
= 


3 2 T 

a 7"2 

dy 


( 1 . 4 . 9 ) 


The conservation of energy equation is integrated across the thickness 
of the film by using Leibniz’s rule and the kinematic condition at the free 
surface. Since the steady- state solution of the fluid velocity will be 
known from the fluid mechanics analysis, the steady solution of the 
conservation of energy equation will be found. The continuity equation 
(1.3.6) and the conservation of energy equation are 


r A V£ = 


(2r) 


^ ^ — constant 


^b _ (2rr) A a 

W = 1T~ 


[91 


Sfl 


by 

6 

3y 

0- 


( 1 . 4 . 10 ) 


Implementing the heat flux boundary conditions corresponding to simple 
heating with no evaporation at the free surface results in the following 
equation. 
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(1.4.11) 


_ (2rr) A q 


W 


Vr 


This equation can be integrated directly to determine the difference in the 
mixed- mean temperature between the inlet and any point in the r- direction. 



1.5 NTTfERICAL SOLUTION PROCEDURE 


(1.4.12) 


Due to the similarity between the hydraulic jump and the shock wave in gas 
dynamics, the MacCormack explicit method (Anderson et al., 1984), which is 
quite often used for the solution of compressible flow problems, will be 
used in the present numerical analysis of this incompressible film flow. 
Since it is an explicit method, the unknown variables are found in terms of 
known quantities, as opposed to implicit methods which must solve a matrix 
equation to obtain, the solution of the problem. 

The governing equation for the present problem is equation (1.3.9): 


v T + = h 


The forward- predictor finite- difference equation for the governing 
equation is forward in time and space. g 


V n+T _ v n 
V k ' V k 


i- 

k H 


Ci - c ; 


+ ArH, 


(1.5.1) 
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The finite- difference equation based on the predicted solution using a 
forward- time, backward- space differencing scheme is: 


{V H7T 


)' 


)X At 

k H 


,n+T r n+i 

J k ‘ G k- 1 


ArH 


The corrected solution is the arithmetic average of the past and 
predicted solutions. 


,n+l 
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(Vf 1 ) 


,n+ 




,n+l 
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P n+i n n+l 
G k ‘ G k-1 


ArH 


nTT 


j (1.5.2) 


Since the forward- predictor velocity is in terms of a forward- space 
approximation, an outlet boundary condition on the velocity is needed. 
For the case of one g, it is assumed that the Froude number at the outlet 
is unity, which is a common boundary condition when a liquid falls over an 
edge because the liquid accelerates from a subcritical velocity to a 
supercritical velocity through the critical velocity. It was 
experimentally observed, however, that for thin films the surface tension 
greatly alters this boundary condition. The Froude number at the inlet and 
the dimensionless velocity at the outlet are related as follows. 


out 


in 


A/3 


Fr. 


_ O 


/3 
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•out 


(1.5.3) 
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For the zero g case and when the Rossby number is Ro < 1, the slope of 
the dimensionless velocity at the last node is set equal to the slope at 
the next to last node for the outlet boundary condition. 

The solution of the governing equation using MacCormack’s method proceeds 
as follows: 


• The parameters pertaining to the numerical domain and the inlet and 
outlet boundary conditions are specified. 

• The initial velocity distribution is input to the program. 

• The variables G and H are computed using the velocity profile at the 
old time step. 

• The velocity distribution at the midpoint time step is calculated in 
terms of the velocity, G, and H at the old time step. An outlet 
boundary condition is needed in this step because of the forward- space 
approximat ion . 

• The variables G and H are computed again by using the velocity profile 
at the midpoint time step. 

• The velocity distribution at the new time step is calculated. The 
inlet boundary condition is used in this step because of the 
backward- space approximation. 


22 


• The values of the velocity distribution at the new time step are used 
as the initial velocity profile for the next iteration. 

• The process is repeated until steady values are reached. 


In explicit schemes, the magnitude of the time step is chosen based 
upon a Courant number, which is defined for the present problem in the 
following manner. 


Ar 

r ( 


V + TU 


= ^ (1 + Ro' 1 ) 


In this study, the grid size was set to a specific value and the time 
step was varied to avoid convergence difficulties. The largest time step 
that did not lead to a divergent solution was then used in the 
calculations. The following table is a general guideline which was used to 
determine the time step size. 


Flow situation 
Plane flow 
Radial flow 

Radial flow with rotation 


Courant number, C 

< 1.0 

< 0.5 

< 0.1 


One hundred grids were used in the plane flow calculations. For radial 
flow, 500 grids were used because the governing equation was very lightly 
damped, so that the numerical oscillations were excessive with only 100 
grids. There was no significant change (< 0.1 percent) in the numerical 
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results when the aforementioned grid specifications were reduced by 
one- half to check the grid independence of the solution. 

1.6 RESULTS AND DISCUSSION 

The thicknesses of free surface films in plane and radial flow have 
been calculated numerically. In addition, the Nusselt number was 
calculated for simple heating with no evaporation at the free surface with 
water as the working fluid. The results of the one- dimensional 
quasi-steady governing equation of fluid motion are given in Figs. (1.3 - 
1.7). The heat transfer results are shown in Fig. 1.8 for plane and radial 
flow in one g and in Fig. 1.9 for radial flow with solid- body rotation. 

Figure 1.3 shows the solution for channel flow as it progresses in 
time. It can be seen that the supercritical portion of the flow develops 
very quickly while the subcritical region and the exact location of -the 
hydraulic jump takes more time to reach the steady state. The Gibb’s 
phenomenon can easily be noticed in the transient solutions of the 
governing equation, but the oscillations become more damped in the steady 
solution.. Due to the assumptions pertaining to the governing equations, the 
transient solutions of the problem are not accurate, but the steady 
solutions should be acceptable. In the following figures, only the steady- 
state solutions are given. It is clear that the velocity profile in the 
immediate vicinity of the hydraulic jump is not a good representation of 
the flow field due to the transient nature of the jump. However, this does 
not affect the numerical results in the supercritical or subcritical 
regions as well as the location of the jump. In order to capture the 
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details of the jump, a two-dimensional model using body- fitting coordinates 
is needed. 

As can be seen in Fig. 1.3, the film thickness increases along the 
plate as the fluid travels away from the inlet of the region. The fluid is 
prevented from rolling back due to the action of gravity by the momentum of 
the fluid in the main direction of the flow. The reason that the film 
thickness increases with distance in the supercritical region is that the 
friction at the interface between the fluid and the plate slows the fluid, 
which translates into an increase in the film thickness due to continuity. 

The effects of changing the inlet Reynolds number and the inlet Froude 
number in plane flow are presented in Fig. 1.4. The graphs m the 
left-hand column of the legend show the dimensionless liquid height when 
the Froude number is held at Fr in - 10.0 and the Reynolds number is varied 
f or Re. n = 25.0, 50.0, and 75.0. By increasing the inlet Reynolds number, 
the location of the hydraulic jump moves downstream and the dimensionless 
height of the liquid at all points in the domain decreases. Physically, 
increasing the Reynolds number can be interpreted as increasing the 
inertial forces on the fluid, so one would intuitively expect that the 
velocity of the fluid would increase as the Reynolds number increases. 
This increase in the fluid velocity is correlated to a decrease in the film 
thickness through the continuity equation. 

The graphs in the right-hand column of the legend in Fig. 1.4 present 

the dimensionless liquid film height when the Reynolds number is constant 

at Re- = 100.0 and the Froude number is varied for Fr- = 5.0, 7.0, and 

in 
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9.0. The location of the hydraulic jump moves downstream when the inlet 
Froude number increases, which has been seen in qualitative experiments. 
The dimensionless height of the liquid in the supercritical region 
generally lies along the same line because the viscous forces are not 
changing in this case. Increasing the Froude number can be thought of as 
decreasing the effect of the gravitational force, so by increasing the 
inlet Froude number, the flow approaches the zero g case. In the zero g 
situation, as shown in Fig. 1.5 for plane and radial flow, the hydraulic 
jump does not occur because the flow at all points in the domain is 
supercritical, since the Froude number approaches infinity in zero g. In 
Fig. 1.5, the dimensionless film height decreases as the Reynolds number 
increases because the inertial forces on the film flow are greater. 

In Fig. 1.6, the inlet Froude number and the inlet Reynolds number in 
radial flow are varied parametrically. Vhen the inlet Froude number is 
held constant and the Reynolds number is increased, the behavior of the 
flow is similar to plane flow in that the film thickness decreases and the 
hydraulic jump moves downstream. The film thickness in the supercritical 
region, however, decreases below the initial thickness downstream from the 
inlet- at high Reynolds numbers. This is due to the increase in the flow 
area downstream as the flow spreads out radially. After this initial 
decrease in the film thickness, the film height increases in the 
supercritical region even though the flow is spreading out radially. This 
is due to the fact that the frictional effect is greater than the effect 
due to the increase in the flow area. The friction at the plate slows down 
the flow and increases the film thickness. As shown by Rahman et al. 
(1989), the liquid film height in the supercritical region may or may not 
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Figure 1.5. Dimensionless liquid height versus dimensionless distance 
r plane and radial flow with f in =0.0 (plane), £ in =10.0 (radial), £ out =100.0 
for different values of the inlet Reynolds number in zero g 





have a minimum value other than the inlet height, depending upon the 
coefficient of friction at the plate and the inlet Froude number. 

The effect of increasing the Froude number while the Reynolds number 
remains constant in radial flow is also shown in Fig. 1.6. The 
supercritical regions of the three cases overlap as in channel flow because 
the viscous forces on the liquid film are not changing. The hydraulic jump 
is moved downstream from the inlet as the Froude number is increased. 

Figure 1.7 presents the effect of solid- body rotation on the 
dimensionless liquid film height in radial flow for one g and zero g. The 
Rossby number of Ro in = 10 10 signifies that the inertial forces are 
dominant in the flow, so that a profile similar to the previous results is 
obtained. For the other cases presented where the Rossby number is Ro- = 
10 and 10 , the hydraulic jump is completely washed out and the 
thickness of the fluid decreases at all points in the domain, which is 
physically realistic. The results for the low Rossby numbers also show 
that the dimensionless liquid film height is the same whether the 
gravitational body force is present or not. This is due to the fact that 
the centrifugal body force in this case is much greater than the 
gravitational body force. 

Figure 1.8 presents the Nusselt number versus the dimensionless 
distance for plane and radial flow under one g. For these calculations, 
the values in equation (1.4.6) were chosen to be = 2.54 x 10 m, 

K = 0.682 V/m-K, v - 2.90 x 10 m^/s, and g = 9.81 m/s^. The values of 
the thermal conductivity and the kinematic viscosity correspond to water at 
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373 K. As would be expected, the heat transfer increases as the Reynolds 
number increases because the height of the film decreases at all points 
along the plate. 

Figure 1.9 presents the Nusselt number for the case of solid- body 
rotation in radial flow under one g. The heat transfer to the liquid in 
this case is two orders of magnitude higher than that without rotation 
because the liquid film becomes extremely thin and the velocity of the film 
is very high. 

It should be noted that while the heat transfer results that were 
presented are for thin films under one g, the heat transfer to thin films 
in a zero g environment can also be calculated with equation (1.4.5). 

1.7 MNCLUSIONS 

In studying plane and radial flows experimentally, it can be seen that 
the flow can experience a hydraulic jump. This is a "mixed" flow situation 
where the flow changes from a supercritical condition to a subcritical 
condition with a sudden increase of the liquid height and a decrease in the 
velocity of the liquid. The similarity between the hydraulic jump and the 
shock wave in gas dynamics suggested the use of the approach where the flow 
is modeled as a transient phenomenon and allowed to march in time to 
achieve the desired steady- state results. Therefore, the equations of 
motion including the transient terms were integrated across the thin liquid 
layer, nondimensionalized, and discretized using the MacCormack explicit 
method to solve for the velocity and height of the liquid film in Cartesian 
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and radial coordinates. Also, the steady- state -heat transfer from the 
* 

plate to the liquid was examined for the case of simple heating with no 
evaporation at the free surface. Since the numerical method of solution is 
explicit, large amounts of computer storage were not necessary. 

For plane flow in one g, a hydraulic jump was predicted which could be 
forced away from the inlet by increasing the initial Froude number or 
Reynolds number, which agrees with what has been found in qualitative 
experiments. In zero g, the liquid film height increased monotonically , so 
that a hydraulic jump was not predicted. 

For radial flow in one g, the liquid film height first decreased for 
high Reynolds numbers due to the increase of the flow area, and then 
increased to form a hydraulic jump. The hydraulic jump in radial flow 
could also be forced away from the inlet by increasing the initial Froude 
number or Reynolds number. In zero g, the hydraulic jump was not predicted 
in radial flow. The effect of solid- body rotation for radial flow in one g 
was found to thin the liquid film and "wash out" the hydraulic jump, which 
is physically realistic. The dimensionless film thickness under one g and 
zero g was the same for solid- body rotation because the centrifugal force 
was much greater than the gravitational body force. The heat transfer from 
the rotating disk was found to increase dramatically when under solid- body 
rotation compared to no rotation because the liquid film is very thin and 
moves with a high velocity. 
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Section II 


EXPERIMENTAL MEASUREMENT OF FILM HEIGHT AND VISUALIZATION OF FLOV 
2.1 SUMMARY 


In this experimental study, the thickness of a thin liquid film with a 
free surface on a stationary and rotating disk is measured with a 
non- obtrusive capacitance technique. The measurements are taken when the 
rotational speed ranges from 0 300 RPM and the flow rate varies from 
7.0 - 15.0 LPM. A photographic study of the thin liquid film was also 
performed and the results are presented herein. When the disk is 
stationary, a hydraulic jump is formed on the disk, which separates the 
supercritical and subcritical regions. As the flow rate of the fluid 
increases, the hydraulic jump is moved away from the center of the disk. 
The liquid film thickness in the supercritical region is affected by the 
inertial and frictional forces on the fluid and by the radial spreading of 
the film across the disk. The film thickness in the subcritical region is 
mainly determined by the surface tension of the fluid because the radius of 
curvature at the outer edge of the disk, which is nearly constant for flow 
rates up to approximately 13.0 LPM. Vhen the disk is rotating, the film 
thickness depends upon the frictional, inertial, and centrifugal forces 
acting on the liquid. In the region near the center of the disk, the 
effects of friction and inertia are greater than that due to the 
centrifugal force, while at the outer edge of the disk the opposite is 
true. A flow visualization study was carried out to examine the nature of 
the free surface of the thin liquid film when the disk is stationary and 
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rotating. When the disk is stationary, surface waves were present in the 
supercritical and subcritical regions at all flow rates studied. The waves 
in the supercritical region have smaller amplitudes than those in the 
subcritical region. At the hydraulic jump, a ’'roller" with a circular 
cross section was found at low flow rates. As the flow rate increased, the 
roller flattened until it merged with the hydraulic jump. The surface 
tension at the edge of the disk held the thickness at this location nearly 
constant, except at higher flow rates where the inertial forces of the 
fluid became greater than the surface tension and decreased the height of 
the fluid at the edge. This effect was also present at low rotational 
speeds, where the surface tension created a standing wave at the edge of 
the disk. As the rotational speed increased, the film changed from the 
wavy- laminar region to a region in which waves ran nearly radially across 
the disk on top of a thin substrate of fluid. These waves appeared to be 
the mechanism through which most of the fluid drained from the disk. 


2.2 INTRODUCTION 


The heat transfer from a rotating disk to a thin film has been 
recognized in the past to be superior to conventional falling films because 
the centrifugal force tends to thin and accelerate the film, which 
translates into higher heat transfer coefficients. As a first step in 
determining the heat transfer to a thin film, it is proposed to 
experimentally study the hydrodynamic characteristics of a thin film with a 
free surface flowing over stationary and rotating disks. The present 
experimental investigation is directed toward the use of a thin film on a 
rotating disk to promote the absorption of a vapor into the liquid on the 
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disk. Specifically, the absorber of a space- based vapor- absorption 
refrigeration system will use a liquid film which is thinned by the 
centrifugal force on a rotating disk to enhance the absorption of the 
refrigerant vapor into the absorbent because a falling film cannot be 
produced in a microgravity environment. The vapor- absorption cycle is more 
appropriate for a microgravity application because the vapor compressor is 
replaced by a liquid pump, which is lighter, requires less maintenance, and 
is easier to manufacture for space- based applications. 

To properly study any problem, the previous studies on the subject 
must be thoroughly investigated. Each of the following researchers have 
experimentally examined the flow of thin liquid films. A brief description 
of their methods and conclusions are provided. 

Vatson (1964) analysed a free- falling jet which impinges on a 
horizontal plane. The jet spreads out radially in a thin him and is 
surrounded by a hydraulic jump, outside of which the thickness of the 
liquid is much greater. Vatson experimentally measured the diameter of the 
hydraulic jump (r^ with dividers, the depth of the fluid outside the jump 
(d) with a point gage, and the flow rate (q) with a measuring jar and stop 
watch. The experimental data ranged from 25.4 < r t < 177.8 mm, 3.30 < d < 
16.5 mm, and 0.73 < Q < 26.8 LPM. The jet Reynolds number Re = Q/av ranged 
from 7 x 10 3 to 1.2 x 10 5 , where a is the radius of the jet. The 
theoretical results of the location of the hydraulic jump and the depth of 
the liquid outside the jump were compared to experimental results with a 
satisfactory agreement. 
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Espig and Hoyle (1965) experimentally jneasured the maximum thickness 
of a liquid film on a horizontal rotating disk. The maximum thickness is 
the height of the film which is measured to the crest of the waves on the 
film. Measurements were taken by using a needle probe which was lowered 
with a vernier, which completed an electrical circuit when the probe 
touched the film surface. The Reynolds number Re = (4 Q/r/iD) ranged from 
10 to 600, where Q is the mass flow rate, n is the dynamic viscosity, and 
D/2 is the distance from the center of the disk. The flow conditions were 
also observed visually with a stroboscope which showed rivulets, 
circumferential waves, and helical waves. The experimental results were in 
agreement with the theoretical and experimental results of previous 
authors. 

Butuzov and Rifert (1972) experimentally measured the heat transfer 
of condensing steam onto one side of a rotating copper disk, the opposite 
side of the disk being cooled with a condenser. The results for the 
average heat transfer coefficient versus the disk angular speed were 
presented. The disk speed varied from 10 to 224 rad/sec and the heat flux 
ranged from 2 x 10^ to 19 x 10^ V/m^. The experimental Reynolds number of 
the condensate was compared to the theoretical results with an agreement to 
within 5-10 percent. The theoretical results were obtained by using a 
previously derived equation for the thickness of a laminar liquid film on a 
rotating disk. This was related to the flow rate of the condensate on the 
disk for steam condensation. 

Charwat et al. (1972) studied the effects of varying the viscosity and 
surface tension of a thin liquid film on a rotating glass plate. The film 
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thickness was determined by the infra-red absorption technique, and it was 
found that the thickness varied from 1 - 150 im for flow rates up to 
0.78 LPI. Concentric, spiral, and irregular waves were found on the 
surface of the film, and were strongly affected by varying the viscosity 
(1 - 2.5 cP) and surface tension (20 - 72 dynes/cm). 

Mat sumo to et al. (1973) compared various theoretical solutions of 
previous authors with a polynomial approximation for the thickness of a 
thin liquid film on a rotating disk. The authors also devised an 
experiment where constant temperature liquid was supplied to the center of 
a rotating disk by a feed nozzle above the disk. The liquid film height 
was measured along the radius by a needle attached to two micrometer 
screws. The kinematic viscosity of the liquid ranged from 9.61 to 58.3 cS, 
the rotational speed varied from 250 to 1500 R.PM, and the flow rate ranged 
from 0.3 to 1.87 LPM. It was concluded that a polynomial of fourth degree 
or higher agreed very well with the experimental results. 

Miyasaka (1974) compared the results of an experimental study with 
those of a theoretical study of the thickness of a thin viscous film on a 
rotating disk. The liquid film was generated with a jet of liquid falling 
onto the center of a horizontal rotating disk. The film height was 
measured by comparing the resistance of the liquid on the disk with that of 
a standard thickness of the liquid, which resulted in an accuracy of * 0.02 
mm. The jet Reynolds number was varied from 200 to 18,000. The theoretical 
values were obtained by solving the governing equations of motion with the 
boundary- layer approximation. The theoretical results include the values 
of the liquid film when viscosity is and is not present. It was concluded 
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that the theoretical values were in agreement with the experimental 
results. 

Ishigai et al. (1977) experimentally measured the liquid film 
thickness and heat transfer from a thin film generated by an impinging 
liquid jet onto a perpendicular surface. The film thickness was measured 
with a needle and micrometer arrangement, and a voltage source applied to 
the surface of the plate. When the needle touched the surface of the film 
an electric circuit was completed, which was read by an oscilloscope. The 
flow rate ranged from 3.0 to 30 LPM. The experimental data of the liquid 
film height was compared to the analytical equation given by Vatson (1964) 
with satisfactory results. 

Labus and DeVitt (1978) experimentally examined the flow patterns of 
the free surface of an impinging jet of liquid on a disk perpendicular to 
the jet in zero gravity. This was achieved in a drop facility in which a 
2.2 second period of 10' 5 g’s could be sustained. Flow visualization 
studies revealed that surface tension and inertia were the major forces 
acting on the liquid. 

Craik et al. (1981) experimentally studied the circular hydraulic jump 
formed by an impinging jet of liquid on a horizontal plate. The liquid 
film thickness in the region near the hydraulic jump was measured using a 
light- absorption technique with a laser and a strong dye. The parameters 
which were varied are as follows: the flow rate, 0.27 to 1.56 LPM; the jet 
radius, 1.0 to 2.15 mm; the jump radius, 12.0 to 40.0 mm; and the outer 
depth, 1.8 to 3.5 mm. The experimental data was compared to Vatson ’s 
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(1964) theory with unsatisfactory results. Flow visualization studies were 
performed, and an eddy just downstream of the hydraulic jump was found. 

Muzhilko et al. (1983) experimentally measured the liquid film 
thickness on a horizontal rotating disk where the liquid is introduced onto 
the center of the disk by an impinging jet. The film thickness was 
measured by sensing electrodes embedded into the surface of the disk, whose 
resistance is a function of the film thickness. The mean film thickness 
was measured at radii of 30, 60, and 90 mm. The angular velocity ranged 
from 95 to 1900 RPM, and the flow rate varied from 0.04 to 1.2 LPM. The 
mean film thickness results were correlated by an empirical equation which 
predicts the data to within ±15 percent in the laminar- wavy region. 

Carper et al. (1986) experimentally studied the heat transfer from one 
side of a rotating disk with an approximately uniform surface temperature. 
The liquid was supplied to the center of the disk by a nozzle. The average 
Nusselt number was presented for the following parameters: the rotational 
Reynolds number (16,000 < Re r = »D 2 n < 545,000), the Prandtl number 
(87 < Pr < 400), and the impinging jet Reynolds number (180 < Rej = dUj/v < 
1300). D is the disk diameter, u is the angular velocity, v is the 
kinematic viscosity, d is the jet nozzle diameter, and Dj is the average 
jet velocity. An empirical equation for the average Nusselt number is 
derived by a multiple linear regression analysis which represents 95 
percent of the data to within ±30 percent. 

Most of the previous experimental investigations were concerned with 
thin liquid films which were generated by an impinging jet onto the center 
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of the rotating or stationary disk. In the present study, the liquid film 
is generated on the rotating disk by a circular collar which directs the 
pressurized liquid onto the center of the rotating disk with a constant 
initial height and velocity. This collar arrangement can be extended to 
collars on each side of the disk, and then multiple disks can be stacked 
onto a rotating pipe so that a very large surface area for absorption can 
exist in a small volume. Also, this arrangement is amenable to numerical 
simulation because of the known inlet conditions. 

The effect of two parameters on the thickness of the liquid film were 
studied: the flow rate (7.0-15.0 LPM) and the rotational speed (0-300 RPM) 
at the steady state. The liquid film thickness was measured with a 
non- obtrusive capacitance probe and digital micrometer arrangement. The 
characteristics of the waves on the free surface were studied with a video 
camera at low and high shutter speeds. 

2.3 EXPERIMENTAL SETPP 


The purpose of this experiment was to measure the liquid film 
thickness on a stationary and rotating horizontal disk where the liquid 
emanates from a pressurized vessel in the center of the disk. In the 
future, this apparatus can also be used to determine the heat transfer from 
the heated disk to the thin liquid film when the disk is stationary and 
rotating. De- ionized water was used as the working fluid because the 
properties are well-known and it is non- toxic. 
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The schematic of the experimental setup j.s shown in Fig. 2.1. The 

f 

specifications of all of the equipment used in this experiment are listed 
in Appendix B following this report. The 406.4 mm dia. rotating disk [1] 
is mounted on a high-precision stainless steel spindle [3] which is hollow 
from the 101.6 ram dia. stainless steel collar [2] to the rotating union 

[21]. The rotating disk is made of aluminum with a surface finish of 

3.8 x 10" 4 mm. The working fluid flows from the circulating pump [25] 
through the stationary pipe to the rotating union, which couples the 
stationary pipe and the hollow spindle. The liquid then passes up above 
the rotating disk and through eight. 3.17 mm dia. radial holes in the hollow 
rotating spindle, which are covered by the collar. The collar directs the 
eight liquid streams down and outwardly so that the fluid is directed onto 
the top side of the disk with a velocity only in the radial direction and 
with a uniform initial height. The collar also makes sure that the height 
of the liquid at the collar does not change with respect to the azimuthal 

direction. The gap height between the bottom face of the collar and the 

top of the disk (0.267 mm) is set by a stainless steel shim which 
separates them. After the fluid flows over the disk, it is collected by 
the water tank [20] to be recirculated. The flow rate of the liquid can be 
varied by the large and small metering valves [26 and 27] , which act as the 
gross and fine adjustments. The volumetric flow rate was determined with 
the turbine flow sensor [28] and the digital flow meter [29] . The 
rotational speed was varied by the frequency inverter speed controller [5] 
attached to the electric motor [4] and is read with the rotary encoder [7] 
and the speedometer [8] . 

In order to decrease the vibrations in the system, it is mounted onto 
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a machine table. The total mass of the frame ^ and the machine table is 
approximately 450 kg, which helps to damp out oscillations. The table is 
also fitted with four vibration mounts which are made of neoprene and are 
approximately 25 mm thick x 75 mm in diameter. The system is located in 
the corner of the building, which aids in decreasing the vibrations seen in 
the laboratory by being very close to the foundation. 

To measure the height of the liquid film, a non- obtrusive measurement 
technique was devised as shown in Fig. 2.1. A non- contact capacitance 
sensor [9] is used to locate the surface of the disk and the surface of the 
liquid. The capacitance sensor directly relates the strength of the 
electric field to the air gap between the sensor and the target. This 
relation is linear within the stated range of the probe. The sensor is 
attached to a digital positioner [11] by which the sensor can be raised or 
lowered. The probe and digital positioner are both mounted onto a linear 
slide [13] which moves the sensor along the radius of the disk. The radial 
location of the capacitance probe is measured by a linear potentiometer 
attached to the linear slide. The procedure to measure the liquid film 
height is described below. The procedure is shown graphically in Fig. 2.2. 

1. The capacitance probe was calibrated to the digital micrometer by 
bringing the probe down until it touched the disk very lightly, which 
could be seen with a light directed behind the probe. This is the 
datum 1 line shown in Fig. 2.2. The capacitance sensor was then 
zeroed at this point with the offset adjustment. The capacitance 
probe was raised slightly above the disk, where the digital micrometer 
was zeroed. The probe was then raised until it read the maximum value 
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o(r) = h +Ab(r) + h ' (r) - h"(r) 


Figure 2.2. Liquid film height measurement 
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of its range. The value read on the capacitance sensor readout [10] 
was compared to what was shown on the digital micrometer readout [12] . 
Any discrepancy could be eliminated by a gain adjustment which is 
supplied on the capacitance sensor readout. This set the slope of the 
linear curve which relates the electric strength to the height of the 
air gap. It was found that the digital micrometer does have a small 
amount of backlash (< 0.008 mm), but this could eliminated by raising 
the probe to the desired location instead of lowering it. 

2. With the disk dry, the capacitance probe was then moved to the datum 2 

line with the digital micrometer, which was set to zero at this 

location. The height at the inner radius was set to an arbitrary 

value h, which is within the range of the probe. 

3. The probe was moved along the radius to measure the deviation of the 

disk from the datum 1, which is called the tare data (h + Ah). 

4. The probe was then raised to datum 3. The liquid flow was started and 

allowed to reach the steady state. The distance between datum 2 and 
datum 3. (h’) was shown on the digital micrometer. The distance from 
datum 3 to the liquid film surface (h M ) was measured by the 

capacitance probe. The liquid film height 6 along the radius of the 

disk can be found with the following equation: 

6 = h’ + (h + Ah) - h" 

The capacitance sensor readout [10] is equipped with a 10 VDC output 
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which is proportional to the air gap between the probe and the target. 

This voltage signal was sampled by a Fluke 2280A Datalogger at a rate of 
one every two seconds to average the readings and to determine the standard 
deviation of each reading. This procedure was necessary due to the waves 
present on the surface of the liquid film. Thus, the liquid film height 
reported here is actually the mean film height. Also, the capacitance 
probe has a sensing spot which has a diameter of 11.28 mm, so a given 
reading was assigned to the radial location underneath the center of the 
sensing spot. Vhen the disk was rotating, the same procedure was followed 
except that the tare data (h ♦ Ah) and the air gap measurements between the 
probe and the free surface (h«) were averaged in the azimuthal direction as 

well. 


The heat transfer from the heated disk to the liquid will be studied 
in future experiments by using a 6 KV etched foil heater [15] between the 
underside of the disk [1] and the ceramic disk insulator [16] . The heater 
is held between the plate and the insulation by screws that are through the 
insulation and are threaded into the plate. The heater is supplied with a 
piece of backing paper that will be placed between the heater and the 
insulation to assure that the heater is evenly compressed against the 
plate. An etched foil heater was chosen for this purpose due to the 
improved uniformity of the heat flux compared to a conventional wire 

heater. 

Since the disk is able to rotate, a sealed slip- ring capsule [19] is 
employed to pass the power from the stationary source [23 and 24] to the 
disk, and to pass thermocouple readings from the disk to the digital 


50 


thermometer [18]. To avoid problems with creating, extra junctions on the 
thermocouples mounted on the rotating disk, the slip- ring capsule was 
manufactured with six constantan wire circuits for thermocouples. In this 
way, an extraneous junction will not be created because a constantan wire 
and a copper wire coming from two separate slip rings will be used to 
create the thermocouple bead. The copper and constantan wires from the 
stationary part of the slip rings will be directly attached to the digital 
thermometer. The only extraneous junction is on the brushes and the gold 
slip rings. Since the brushes and the gold slip rings are very small, these 
junctions are nearly isothermal, so_ the contribution due to these junctions 
should be negligible. Thermocouples [17] mounted in holes in the disk 
along one radius will monitor the temperature of the disk while other 
thermocouples will measure the temperature of the liquid before and after 
it flows over the disk. The thermocouples in the holes in the plate will 
be a distance of 1.58 mm from the surface of the disk in order to measure 
the interface temperature as closely as possible. The thermocouples will 
be cemented into the holes with a thermally conductive epoxy which will not 
be degraded by the operating fluid. The mean temperature of the fluid at 
different radial locations will be measured by a sheathed thermocouple 
which will be moved with the digital micrometer [11] and the linear slide 
[13]. Since the liquid is in a closed system, a heat exchanger [22] will 
be used to keep the temperature of the inlet water at a constant value. 

2.4 ERROR, ANALYSIS 

The errors in the liquid film measurements are described and quantified, 
which involve two instruments: the digital micrometer, and the capacitance 
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sensor. A brief summary of the method to determine the error will now be 
given. A more detailed description is presented in Appendix C and a full 
discussion on error analysis is given by Miller (1989). 

The liquid film height is given by the following equation shown in 
Fig. 2.2: 

6 = h’ + (h + Ah) - h" (2.4.1) 


Following the standard error analysis procedures, the root- sum- square error 
of the system is given by 


E 


rss 


(Ah’)2 + [A(h + Ah)]^ + (Ah M )^ 


(2.4.2) 


The first term in equation (2.4.2) is the stated accuracy of the 
digital micrometer, which was calibrated at the factory with an NBS 
traceable Mark- Tech Laser Gage Model 7980. The maximum deviation over the 
range of movement reported on the calibration certificate is 0.001 mm, 
which was taken as the error of the instrument. 


The second term involves the accuracy of the capacitance sensor, which 
is given by the manufacturer to be ±0.1 percent of the range when 
calibrated to a known standard, which was the digital micrometer. The 
capacitance probe was calibrated at the beginning of each test against the 
digital micrometer as described in Sec. 2.3. The maximum deviation from 
linearity of the capacitance probe was 0.008 mm. This value was added to 
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the stated accuracy of the probe. 

The second term in equation (2.4.2) takes on different values when the 
disk is stationary or rotating. Vhen stationary, the second term is the 
accuracy of the capacitance sensor given by the manufacturer plus the 
maximum deviation from linearity, which is 0.01 mm. 

Vhen the disk is rotating, the error due to the disk not being 
perfectly flat and the wobble caused by the bearings must be taken into 
account. This error is given by the sum of the standard deviation of the 
air gap measurement between the probe and the dry disk taken at the outer 
edge of the disk (where it is maximum) and the error of the capacitance 
sensor. Multiple populations with an increasing number of samples were 
taken and the means were compared with a two-tailed z-test to give a sample 
number which ensures a mean at the 95 percent confidence level. The 
standard deviation of the dry rotating disk with 50 samples is a - 0.02 mm, 
so the error in the second term in equation (2.4.2) is 0.03 mm. 

The last term in equation (2.4.2) is again the error associated with 
the capacitance sensor. Due to the fact that the mean film thickness is 
being reported, further comments concerning the third term in equation 
(2.4.2) are necessary. At all flow rates and rotational speeds, waves 
were present on the free surface of the liquid film. A test was performed 
to compare the means of different numbers of samples in the subcritical and 
supercritical regions to determine a suitable number of data points for a 
mean which falls in the 95 percent confidence level. For subcritical flow, 
a normal z-test was performed, and it was found that for 125 and 250 
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samples, the means were equal at the 95 percent confidence level, so 250 
samples were taken for each data point reported. Likewise, for 

supercritical flow 80 samples were taken in order to achieve the same 
confidence level. A similar test performed when the disk was rotating 
showed that the number of samples needed for the 95 percent confidence 
level was the same as that for the stationary supercritical region. The 
root* sum- square error in the liquid film height measurement is as follows: 

Stationary disk - E = ±0.01 mm 

rss 

Rotating disk - E rss = ±0.03 mm 

The gap between the collar and the disk was measured with feeler gages 
and the error in this data point was found to be ±0.01 mm. 

2.5 RESULTS AND DISCUSSION 

2.5.1 Liquid Film Thickness Measurements 

The mean thickness of a liquid film of de- ionized water as it flows 

across a horizontal disk has been measured experimentally. The flow rate 

ranged from 7.0 - 15.0 LPM and the rotational rate varied from 0 - 300 

RPM. - Figures 2.3 - 2.8 present the film thickness when the rotational 

speed is held constant and the flow rate is varied. Figures 2.9 - 2.11 

show the film thickness when the flow rate is held constant and the 
rotational speed is changed. In all of the measurements presented in this 
report, a duplicate measurement was made and the two runs were within 5 
percent of each other, which is within the experimental error of the film 

thickness measurement system. 
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Figure 2.3 shows the film thickness when the disk is stationary in 
both the supercritical and subcritical regions. As the flow rate 
increases, the hydraulic jump is moved toward the outer edge of the disk. 
The hydraulic jump is the sudden increase in the film thickness and is 
actually a momentum balance between the supercritical flow and the 
subcritical flow. The exact shape of the jump could not be determined with 
the present measuring instrument because if an attempt was made to measure 
the film thickness in the supercritical region very near the jump, the free 
surface of the subcritical region would touch the side of the capacitance 
sensor and saturate the signal. Therefore, the straight line without data 
points between the supercritical and subcritical regions is where the 
hydraulic jump resided. 

The first data point at 50.8 mm is the gap height between the collar 
and the rotating disk (0.267 mm). The next data point at 76.2 mm is the 
first point which could be read by the capacitance probe. Therefore, the 
line connecting these two points is simply a linear interpolation and no 
conclusions can be drawn concerning this region. Even though the outer 
radius of the disk is 203.2 mm, measurements are not reported past 195.6 mm 
because the sensing spot of the capacitance sensor was over the edge of the 
disk. 


In Fig. 2.3, it can be seen that the film is nearly horizontal in part 
of the subcritical region when the flow rate is 7.0 LPM and then drops off 
near the edge of the disk. This decrease in the film thickness at the edge 
is due to the acceleration of the flow as the liquid exits the disk. The 
radius of curvature at the edge was nearly constant up to approximately 
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Radial location (mm) 

Figure 2.3. Film thickness versus radial location 
for different flow rates (stationary disk) 






13.0 LPM and produces significant surface tension effects. This can be 
noticed in Fig.' 2.3 by the curves showing 7.0 - 13.0 LPM, which approach 
the same film thickness near the edge of the disk. For 15.0 LPM, the 
hydraulic jump was nearly pushed off the edge of the disk, so it did not 
reach the same height as the curves with lower flow rates. Figure 2.3 also 
shows that the liquid film thickness in the subcritical region is an order 
greater than that in the supercritical region. 

In Fig. 2.4, the film thickness in the supercritical region when the 
disk was stationary can be seen to increase along the radius of the disk. 
In the supercritical region, the major forces on the liquid film are those 
due to inertia and friction. Since the frictional forces tend to slow the 
liquid, the film thickness increases due to continuity. 

Figure 2.4 also shows that the film thickness in the supercritical 
region may increase raonotonically or may first decrease and then increase 
downstream. This decrease in the film thickness is due to the radial 
spreading of the fluid as it travels across the disk. This effect is then 
overtaken by that due to friction, at which point the film thickness 
increases. For flow rates of 11.0 LPM or less, the effect of friction is 
dominant over the effect of the spreading of the liquid, so the liquid 
increases monotonically . As the flow rate (or Reynolds number) increases, 
the film thickness decreases because the effect of the radial spreading of 
the film becomes greater than the effect of friction. 


The film thickness for different flow rates when the rotational speed 
is 55 RPM is shown in Fig. 2.5. Up to approximately 120 mm the flow has 
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Radial location (mm) 

Figure 2 . 4 . Film thickness versus radial location in the supercritical 
region for different flow rates (stationary disk) 




59 


Radial location (mm) 

Figure 2.5. Film thickness versus radial location 
for different flow rates (55±2 RPM) 


as the 


the same trend as the supercritical flow of the stationary disk: 
flow rate increases, the film thickness decreases due to the effect of the 
radial spreading of the flow overtaking the effect of friction. Past this 
region, however, the profile takes on a very different appearance. For 7.0 
LPM, the film thickness begins to flatten out and then decrease past 120 
mm. When the disk is rotating, the dominant forces on the film are as 
follows: frictional, inertial, and centrifugal. Near the center of the 

disk, the frictional and inertial forces are prevalent. Closer to the edge 
of the disk the centrifugal force comes into play. Therefore, near the 
edge of the disk one would expect that the film would accelerate because of 
the centrifugal force and become thinner due to continuity. It can also be 
seen in Fig. 2.5 that as the flow rate increases, the point where the 
curves begin to flatten due to centrifugal force travels downstream. This 
is because the effect of inertia on the flow is dominant over the effect of 
the centrifugal force for a longer distance. Near the edge of the disk, 
the film thickness increases dramatically due to the effect of surface 
tension at the edge. This is actually a fractional hydraulic jump which 
does not reach its full height before the edge of the disk. As the flow 
rate increases, the height of this jump decreases because the force due to 
the inertia of the liquid overcomes the force due to surface tension. 


Figure 2.6 presents the film thickness measurements for different flow 
rates when the rotational speed is 100 RPM. In this graph, it can be seen 
that the hydraulic jump near the edge of the disk is not present because 
the centrifugal force is greater than the force due to surface tension. 
Also, the curves have the same trend as in Fig. 2.5 where the film 
thickness first increases because of friction and then decreases due to the 
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Figure 2.6. Film thickness versus radial location 
for different flow rates (100±2 RPM) 


centrifugal force. A further examination of jthis figure shows that the 
* 

curve of 15.0 LPI lies above those of 11.0 LPM and 13.0 LPI in the region 
where the film thickness is increasing. A possible reason for this is that 
the difference between these curves is on the order of the experimental 
error in the measurement system. 

Figures 2.7 and 2.8 present the film thickness when the flow rate is 
varied from 7.0 to 15.0 LPI for rotational speeds of 200 and 300 RPM, 
respectively. Again, the same trend is evident where the film thickness 
first increases due to frictional effects and then decreases due to the 
effect of centrifugal force. In comparing these two figures, it can be 
seen that the radial location where the effect of centrifugal force begins 
to dominate the flow moves toward the center of the disk as the rotational 
speed increases. This is shown further in Figs. 2.9 - 2.11 where the flow 
rate is held constant (7.0, 11.0, and 15.0 LPM, respectively) and the 
rotational speed is varied from 55 to 300 RPM. In Fig. 2.9, as the 
rotational speed increases the maximum thickness occurs closer to the 
center of the disk. This can also be seen in Figs. 2.10 and 2.11. Also 
shown in Figs. 2.9 - 2.11 are the supercritical regions for the different 
flow rates when the disk is stationary. It can be seen that the film 
thickness in the supercritical region is generally close to the case when 
the rotational speed is 55 RPM. The exception is near the end of the 
supercritical region where the film thickness increases, as shown in Figs. 
2.10 and 2.11. This phenomenon is due to the fact that the stationary film 
is not being acted upon by centrifugal forces, so the film thickness 
continues to increase in this region because of the effect of friction. 
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Figure 2.8. Film thickness versus radial location 
for different flow rates (300±2 RPM) 
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Figure 2.9. Film thickness versus radial location 
for different rotational speeds (7.0±0.3 LPM) 
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Figure 2.10. Film thickness versus radial location 
for different rotational speeds (ll.0±0.3 LPM) 
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Figure 2.11. Film thickness versus radial location 
for different rotational speeds (15.0±0.3 LPM) 




2.5.2 Flow Visualization 


A photographic study was performed to determine the basic phenomena 
which occur when the disk was stationary and rotating. Vhen the disk is 
stationary, waves can be seen both on the subcritical and supercritical 
regions. The waves on the subcritical region have larger amplitudes than 
those in the supercritical region. At low flow rates, a toroidal "roller" 
was found at the hydraulic jump when the disk is stationary. At higher 
flow rates, the roller flattened out until it could not be distinguished 
from the increase in the film thickness at the hydraulic jump. Also for 
high flow rates, the surface tension of the liquid at the free surface 
creates a standing wave at the edge of the stationary disk. This standing 
wave is also present at all flow rates when the disk is rotated at low 
angular speeds. At higher spin rates, waves appear on the disk which flow 
nearly radially across the disk on top of a thin substrate of fluid. 

In Figs. (2.12 and 2.13), the disk was stationary, the flow rate was 
set to 7.0 LPM and the shutter speed was 1/1000 second. Vith this shutter 
speed, the waves on the free surface of the liquid film could be 
photographed.. In all of the photographs presented, the direction of the 
flow is from right to left and the direction of rotation is clockwise. The 
waves in the supercritical region are shown in Fig. 2.12. Near the collar 
the amplitude of the waves is small, but close to the hydraulic jump the 
amplitude of the waves becomes larger. The waves in the subcritical region 
can be seen in Fig. 2.13. These waves have a much larger amplitude than 
those in the supercritical region. This was verified when the liquid film 
height was measured with the capacitance sensor. Even though the actual 
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height of the waves could not be measured, .the standard deviations of the 

t 

signals from the mean film thickness were much larger in the subcritical 
region than those in the supercritical region. 

Figures (2.14, 2.15 and 2.17) show the behavior of the hydraulic jump 
for low flow rates. These photographs were taken with a shutter speed of 
1/60 second. Figure 2.14 shows the hydraulic jump at a flow rate of 3.0 
LPM. The hydraulic jump occurs immediately after the liquid exits the gap 
between the collar and the disk. The flow rate was increased to 4.0 LPH, 
which is shown in Fig. 2.15. The hydraulic jump is pushed away from the 
collar approximately 30 mm. The shape of the jump can be seen to be a 
toroidal "roller". The motion of the roller is very chaotic in that its 
width and height oscillate in a seemingly random fashion. At this flow 
rate, however, the cross section of the roller is nearly circular. At 
higher flow rates, it was found that the shape of the roller changes 
significantly. A more detailed visual study was performed, and it was 
found that at flow rates between 3.0 - 6.0 LPM, the roller had a nearly 
circular cross section as shown in Fig. 2.16a. For flow rates between 7.0 
8.0 LPM, the cross section of the roller was oval in shape as shown in 
Fig. 2.16b. . For flow rates at and above 9.0 LPM, the shape of the roller 
was no longer distinguishable from the gradual increase in the film height, 
which is presented in Fig. 2.16c. The surface of the film at the hydraulic 
jump at this flow rate could clearly be seen to be falling down toward the 
center of the disk in the direction opposite to the main direction of the 
flow. Therefore, the roller had submerged such that it could not be seen 
photographically, but its presence was unmistakable. This change in the 
shape of the roller is thought to be due to the change in the balance 
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Roller 



Figure 2.16a. Roller for 3. 0-6.0 LPM 





Figure 2.16b. Roller for 7.0— 8.0 LPM 



Figure 2.16c. Roller for 9.0-15.0 LPM 
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between the surface tension and the momentum of the fluid. At lower flow 
rates, the velocity of the fluid within the roller and the surface velocity 
of the liquid film do not have sufficient momentum to overcome the surface 
tension which holds the circular shape of the roller. As the flow rate and 
the surface velocity of the supercritical liquid film increases, the effect 
of momentum of the liquid overcomes that of surface tension, which results 
in the change in the shape of the roller. 

Figure 2.17 shows the film surface when the flow rate is 5.0 LPM. A 
series of small concentric waves immediately upstream of the hydraulic jump 
can be seen, which are illustrated in Fig. 2.16. These concentric waves 
were present at flow rates between 4.0 - 15.0 LPM and had larger amplitudes 
at higher flow rates. 

Figures (2.18 - 2.24) present the rotating disk for a flow rate of 7.0 
LPM and rotational speeds of 55, 100, 150, 200, 250, 300, and 350 RPM, 
respectively. Figure 2.18 shows the disk rotating with a speed of 55 RPM. 
The waves on the free surface appear to be similar to the waves found in 
the subcritical region on the stationary disk with the same flow rate (see 
Fig. 2.13). Measurements of the standard deviation from the mean film 
thickness, however, showed that the amplitude of the waves were closer in 
height to those in the supercritical region. Near the edge of the disk is a 
standing wave created by the surface tension of the fluid. In this case, 
the rotational speed is slow enough that the fluid follows the edge of the 
disk downward instead of flying off the disk horizontally. Therefore, the 
radius of curvature of the free surface around the edge of the disk acts to 
increase the film height immediately upstream from the edge of the disk. 


73 



ORIGINAL PAGE IS 
OF POOR QUALITY 





ORIGINAL PAGE IS 
OF POOR QUALITY 




75 




ORIGINAL PAGE IS 
OF POOR QUALITY 



Figure 2.21. Rotating disk, 7.0 LPM, 200 RPM 
Figure 2.22. Rotating disk, 7.0 LPM, 250 RPM 
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Figure 2.23. Rotating disk, 7.0 LPM, 300 RPM 
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Two flow regimes were found to exist on ,the rotating disk: 

» 

wavy- laminar and radial wave flow. The transition from the wavy- laminar 
flow to the radial- wave flow has been visually examined for 7.0 and 13.0 
LPM and is presented in Figs. (2.18 - 2.31). The wavy- laminar flow is seen 
at lower rotatinal speeds and the radial- wave flow occurs at high 
rotational speeds. The wavy- laminar flow is shown in Fig. 2.18 for 7.0 LPM 
and 55 RPI. In this regime, no recurring wave patterns are distinguishable 
except the end effects at the outer edge of the disk. The radial- wave 
regime is shown in Fig. 2.30 for 13.0 LPM and 300 &PM. Several 
well-defined radial waves can be sepn which appear to flow across the disk 
on top of a thin substrate and carry the bulk of the fluid off of the disk. 
In Fig. 2.30, a radial wave is exiting the disk on the left-hand side in 
the center. Immediately below this wave, very little fluid is leaving the 
disk where the thin substrate is present. The transition between these two 
regimes will now be discussed. 

When the rotational rate was slowly changed from 55 RPM to 100 RPM for 
7.0 LPM, the waves on the free surface began to form a pattern which could 
be easily distinguished, as shown in Fig. 2.19. Parts of the free surface 
broke, free from the wavy- laminar regime to form V-shaped waves that ran at 
a diagonal angle between the radial and circumferential directions opposite 
to the rotation of the disk. At the lower rotational speeds, these waves 
were present only at the outer edge of the disk. The film thickness at the 
base of the ’V’ was greatest, and that at the top of the V was least. It 
is felt that the effect of the centrifugal force eventually overcame that 
of the surface tension, so that the liquid film began to run in rivulets 
toward the edge of the disk. As the rotational speed increased, the width 
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of the V-shaped waves decreased and the length increased until the waves 
traversed the radius of th disk as shown in Fig. 2.21. It can be seen in 
Fig. 2.22 that the spacing between the radial waves increases with the 
radius because of the increase in the flow area. For rotational speeds at 
and above 200 RPM, the radial waves could be seen in almost every still 
photograph, so the transition region for 7.0 LPM is 100-200 RPM. 

Figures (2.25 - 2.31) present the surface waves for 13.0 LPM at 
rotational speeds of 55, 100, 150, 200, 250, 300, and 350 RPM. Dpon 
examining these figures, the transition region was determined to be 
approximately 150 - 200 RPM, which is higher than that of 7.0 LPM. It is 
thought that at the higher flow rate, the transition is delayed because the 
effect of centrifugal force must overcome the increased effect of the 
radial momentum of the fluid. 

Figures (2.18 and 2.25) compare two flow rates (i.e., 7.0 and 13.0 
LPM, respectively) at the same rotational speed (i.e., 55 RPM). In Fig. 
2.18, the width of the standing wave at the edge is much larger than that 
of Fig. 2.25, which has a higher flow rate. This shows that at higher flow 
rates, the momentum of the flow overcomes the effect of surface tension. 

2.6 CONCLUSIONS 

The characteristics of a thin liquid film with a free surface on a 
stationary and rotating disk have been examined experimentally. The film 
thickness was measured for different flow rates and rotational speeds with 
a non- contact capacitance technique. Also, the nature of the waves on the 
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Figure 2.25. Rotating disk, 13.0 LPM, 55 RPM 
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Figure 2.27. Rotating disk, 13.0 LPM, 150 RPM 


Figure 2.28. Rotating disk, 13.0 LPM, 200 RPM 
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Figure 2.29. Rotating disk, 13.0 LPM, 250 RPM 


Figure 2.30. Rotating disk, 13.0 LPM, 300 RPM 
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Figure 2.31. Rotating disk, 13.0 LPM, 350 RPM 
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free surface was determined photographically.^ Vhen the horizontal disk was 
stationary, the thin film experienced a hydraulic jump, in which the 
velocity decreases and the thickness increases at a certain radial 
distance. In the region upstream of the jump, the film thickness was 
governed by the inertial and frictional forces on the fluid, as well as the 
radial spreading of the fluid. Downstream from the jump, the film 
thickness was mainly determined by the radius of curvature of the liquid at 
the outer edge of the disk. The liquid film thickness on the rotating disk 
was affected by the inertial and frictional forces on the fluid near the 
center of the disk, and by the centrifugal force near the outer edge of the 
disk. The flow visualization study revealed the presence of a "roller" at 
the hydraulic jump on the stationary disk, whose shape changed with the 
flow rate. Also, the transition between wavy- laminar flow to radial- wave 
flow was observed on the rotating disk as the rotational speed was 
increased. The range of transition shifted to higher rotational speeds as 
the flow rate increased. 
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Section III 


TVO- DIMENSIONAL COMPDTATION OF THE FREE SURFACE FLOV 
OF A THIN LIQUID FILM USING "PRESSURE OPTIMIZATION METHOD" 

3.1 SUMMARY 

The results of numerical computations are presented for the free 
surface flow of a thin liquid film in the presence or absence of a 
gravitational body force. Three different flow systems were studied: (a) 

a falling film down a vertical wall, (b) plane and radial film flows under 
zero gravity and (c) plane and radial film flows along a horizontal plate 
in the presence of gravity. In the case of film flow along a horizontal 
plate where gravity acts across the thickness of the film, the Froude 
number, which characterizes the flow regime (i.e., supercritical or 
subcritical) , is found to be the most dominant parameter. The 

transformation of the flow from supercritical to subcritical is associated 
with a hydraulic jump. The distributions of the film height, film velocity 
and friction coefficient are presented. 
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3.2 INTRODUCTION 


The flow of a thin liquid film is encountered in many engineering 
devices. A thin film is quite commonly found during evaporation or 
’condensation on a solid surface in a compact heat exchanger or cooling 
tower, spin coating in metal industries, and impingement cooling of a solid 
wall with a liquid jet. Besides practical applications, the fluid 
mechanics of thin film flows is important from a theoretical point of view 
since both viscosity and free- surface effects are significant in these 
flows. Moreover, the understanding of such flows under reduced or zero 
gravity is essential for the proper design of heat exchangers and heat 
pumps for space applications, which was the primary motivation for the 
present study. 

The falling of a thin liquid film along a plane vertical wall has been 
studied by many investigators since the turn of this century. For steady 
fully- developed laminar flow, a theoretical solution can be derived from a 
simple balance between the gravitational body force and the shear force at 
the solid wall (Bird et al. (I960)). The film height remains constant and 
the velocity profile across the film becomes parabolic in the fully 
developed region. The results of developing flow when a film is introduced 
at its equilibrium height is presented in the review aritcle by Faghri and 
Payvar (1979). This review also included the experimental studies on 
laminar flow with constant thickness. A film falling under the influence 
of gravity ceases to be laminar and constant in thickness when the flow 
rate is high. Waves tend to appear on the surface and the flow becomes 
turbulent as the flow rate is increased. A number of theoretical as well 
as experimental studies have been performed to understand the wavy- laminar 
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and turbulent flows (see, for example, Seban and Faghri (1976, 1978) and 
Hirshburg and Florschuetz (1982)). 

A somewhat less studied problem is the spread of a liquid film over a 
plate. Vatson (1964) presented results of analytical and experimental 
studies of the radial spread of a liquid jet impinging on a horizontal 
plane for laminar and turbulent flows. By using the boundary layer 
approximations for the governing equations, analytical solutions using a 
similarity transformation along with the Pohlhausen integral method were 
derived. The analysis covered the regions where the boundary layer 
thickness is less than the film height and where the film is totally 
engulfed by the boundary layer. The effects of the gravitational pressure 
gradient was discussed. The possibility of a hydraulic jump in such a flow 
was also anticipated. However, the analysis was applicable only to the 
supercritical flow before the jump. An equation was also presented to 
predict the jump height for any given location of the jump. The agreement 
between the experimental data and the analysis was satisfactory. 

Another interesting problem of thin film research is the spreading of 
the film under the action of centrifugal force as seen in a rotating 
system. An approximate analytical solution for laminar flow on a rotating 
disk was developed by Rauscher et al. (1973). An asymptotic expansion 
technique was used where the radial spread of the fluid was perturbed to 
determine the effects of convection, Coriolis acceleration, radial 
diffusion, surface curvature and surface tension. These higher order 
effects were discussed on a physical basis. 
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In all of the previous studies concerning thin liquid films, the 

A* *' 

investigators have tried to develop analytical models or have taken 
experimental data. Some of these models are quite approximate in nature 
and do not bring out the finer details of the flow field. Moreover, in all 
of the previous numerical studies (Faghri and Payvar (1979)) concerning the 
laminar falling film, a constant thickness was used in the calculation 
domain and no degree of freedom was permitted for the variation of free 
surface along the flow. A general numerical finite- difference solution of 
a thin film flow accounting for the variation of free surface height is not 
available at the present time. These flows are difficult to solve by the 
finite- difference method since the geometry of the free surface is not 
known ahead of time, and the surface profile cannot be fitted in a regular 
Cartesian or cylindrical coordinate system. Moreover, none of the studies 
mentioned above has considered the flow under reduced or zero gravity, 
which is expected to be different from the flow under normal gravity. A 
proper understanding of such flows is essential in the design of space 
cooling systems. 

The present study is undertaken to develop a general numerical 
solution procedure for free surface thin film flows which can be applicable 
to both plane and radial systems, and to both normal and zero gravity 
environments. The results highlight the effects of gravity for different 
configurations of the flow. 

3.3 PROBLEM FORMULATION 

The equations governing the conservation of mass and momentum in a 
thin film of fluid which is Newtonian with constant properties are given by 
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These governing equations have to be supplemented with appropriate 
boundary conditions. At the solid wall, the no- slip condition exists, 

therefore, V = 0. On the free surface, the shear stress vanishes which 

implies dY / dn = 0, where n is the coordinate normal to the free surface. 

Moreover, in the absence of any significant surface tension, the static 

pressure on the free surface must equal the ambient pressure. By setting p 

equal to the difference between the actual and ambient pressures, then p = 

0 on the free surface since pressure is a scalar quantity. Boundary 

conditions must also be assigned in the direction of the flow at two 

locations: the inlet and exit to the control volume or computational 

domain. The appropriate conditions were determined by analyzing the 

characteristic behavior of the flow, which is presented in a later section. 

These depend on whether the flow is supercritical, subcritical, or mixed. 

For a supercritical flow, which includes most of the cases considered here, 
-+ -* 

h = h . and Y = V- at the inlet and <9V/#n = 0 at the exit, where n is the 
coordinate normal to the exit plane. For subcritical flow, h = h , at the 
exit was prescribed instead of the inlet height. Moreover, the pressure 
was prescribed at the exit boundary. The boundary- fitted coordinate system 
used here is shown in Fig. 3.1 and the boundary conditions in component 
form are listed in Table 3.1 for plane and radial flows. The coordinate 
system is dicussed in detail in a later section of the paper. 




Table 3.1: Boundary Conditions for Plane Flow and Radial Flow 


at y = 0: 


v = v = 0 


at y = 6 : 

at z = 0 or r. : 
in 

at z = L or r out : 



I ' 


t = 0 


v = 0, 


v = V 


in’ 


v = 1.5 V in 


for unifora entrance 
2 


[ 2(£) - (£) ] » for parabolic 

6 entrance 


p = 0, for cases 1 and 2 
P = Pg • y), for case 3 
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The flow configurations considered ia. the present investigation are 
shown in Fig. 3.2. They can be broadly classified into three groups 
according to the presence and orientation of the gravitational body force: 

(1) Falling film along a vertical wall 

(2) Film flow under zero gravity 

(3) Film flow along a horizontal plate in the presence of gravity. 

The three possible orientations with respect to gravity will be denoted by 
1, 2, and 3. Moreover, P and R will denote plane and radial flows and A, B 
and C will denote different combination of flow parameters. The parameters 
used for this study are listed in Table 3.2. 

The first case is a classical problem where the major driving 
mechanism is the gravitational body force. This problem was used to check 
the accuracy of the present numerical methodology. The effects of 
introducing the film at a height other than the equilibrium height for a 
given flow rate will also be investigated for this case. 

For film flow under zero gravity, two problems will be considered: 

(2P) Plane film flow under zero gravity 
(2R) Radial film flow under zero gravity 

In the absence of any gravitational body force, the orientation of the 
plate becomes immaterial. The flow remains the same whether the plate is 
vertical, horizontal, or inclined. In this case the flow is driven by 
inertia and viscous forces. In the radially spreading flow, the area 
available for the film increases downstream and acts as an added mechanism 
for the reduction of the film velocity. 
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2P. Plane flow at zero gravity 2R. Radial flow at zero gravity 



3P. Horizontal plane flow 3R. Horizontal radial flow 


7 ig. 3.2. Flow systems in the present investi 
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Table 3.2: Flow Conditions and Free Surface Specifications 


Case 1PA: = <5 fd , (= 0.-000595 m) 

Re = 12.5, g/g =1, L = 0.03 m 
s 

Free surface: A = 1 

Case 1PB: ^ = 1.26^, <^ d = 0.000595 m) 

Re = 12.5, g/g = 1, L = 0.03m 
s 

Free surface: f 4 = ^ + W**’ for £ < ^1 
1 A = 1, for (> 

Case 1PC: $ in = 0.8<5 fd , ($ fd = 0.000595 m) 

Re = 12.5, g/g =1, L = 0.03 m 

Free surface: { 1 » ( 1 * *<>“ 

* A = 1, for ( > ^ 


Case 2P: 6- = 0.000595 m 

in 

Re = 12.5, g/g = 0, L = 0.03m 

u 

Free surface: A = (1 + A£) an 

Case 2RA: = 0.005 m, r^ n = 0.050$ m, r Qu ^. = 0.1953 m 

Re in = 404, g/g s = 0 

Free surface: A = (1 + A<f) an 

Case 2RB 5 ■ = 0.000508 m, r. = 0.0508 m, r. = 0.1 m 

in ’ m ’ out 

Re- = 8.5, g/g 0 = 0 
in ’ blb s 

Free surface: A = (1+ A<f) an 


Continued on next page 
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Case 3P : 


Re = 88.75, g/g = 1, L = 0.1445 m 

n • _ J 5.661, for case 3PA 

in { 8.582, for case 3PB 

Fr out = 1 -° 

Free surface: 

A = (1 + A^)* 11 , supercritical 
A = C [2 - (D + B<f)k n ] ? subcritical (( t 1) 


Case 3R: 


Re in = 50.5, g/g s = 1, r in = 0.0508 m, r Qut = 0.1953 m 

p _ f 7.442, for case 3RA 

in \ 10.96, for case 3RB 
Fr 0„t * !■« 

Free surface: 

A = (1 + A<f) an , supercritical 
A = C [2 - (D + B^)^ n ], subcritical (£ ^ 1) 
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In the third group, we will again consider two. problems according to 

t 

the flow geometry under consideration: 

(3P) Plane film flow along a horizontal plate in the presence of gravity 
(3R) Radial film flow along a horizontal plate in the presence of gravity 

In the case of horizontal thin film flow where the gravitational body 
force acts across the thickness of the film, different flow regimes 
(i.e., supercritical, subcritical, or both) may be present according to the 
local film velocity and height. The transition of the flow from 
supercritical to subcritical can, take place through a hydraulic jump. 
Since a sudden transition of the flow takes place across a jump, special 
analytical and computational tools are required to calculate the flow 
around this singular point. 

3.4 NUMERICAL SOLUTION PROCEDURE 

The complete governing transport equations (3.1 and 3.2) along with 
the appropriate boundary conditions (Table 3.1) were solved numerically 
using a finite- difference scheme. Since the free surface geometry cannot 
be handled very well with a regular rectangular or cylindrical coordinate 
system, a boundary- fitted curvilinear coordinate system had to be used. In 
this system, the free surface of the film was used as one of the boundaries 
of the control volume. 

A curvilinear system can be either orthogonal or non- orthogonal 
depending on whether the faces of the control cells are orthogonal to each 
other or not. The orthogonal system has the advantage of simplicity 
compared to the non- orthogonal system. In either system, the vectorial 
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form of the governing equations [i.e., eqns.^ (3.1), and (3.2)] can be 
written in terms of components and can be discretized to determine the 
finite- difference equations. In most of the computations presented, the 
coordinate system was non- orthogonal. 

Vithin the range of the general non- orthogonal coordinate system there 
exist several options in formulating the equations. These options arise 
from the freedom available in the choice of velocity components and their 
direction with reference to the coordinates. Thus, velocity and force 
vectors can be resolved either into their Cartesian, covariant or 
contravariant components. Moreover, the problem can be solved in a 
physical domain or transformed into a domain where the grid cells are 
rectangular and other physical quantities are non-dimensional or reduced in 
dimension. Although all these options are obviously equivalent to each 
other from the physical point of view, they are substantially different as 
far as numerical treatment is concerned, each presenting its own problems. 
In the present study, the problem was solved in the physical domain where 
covariant velocity components were used. These are components parallel to 
the cell faces. 

The grid system used can be considered as a distorted version of the 
usual orthogonal Cartesian grid system in which grid lines and control 
cells are stretched, bent and twisted in an arbitrary manner, subject to 
the cells retaining their topologically Cartesian character. This means 
that grid cells always had four sides and four corners in the 
two-dimensional domain considered here. 
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As shown in Fig. 3.1, the local coordinates aje defined along lines 

i 

joining adjacent cell centers. The z-axis was taken in the streamline 
direction and y-axis in the direction across the film. The resolutes of 
the velocity vectors in the y and z directions are v and w, respectively, 
and can be defined as 


v 



j 


w 



• k 


Here j and k are unit vectors in the direction of the coordinate axes. In 
general, the resolutes are not the same as the velocity components in these 
directions, but can be related to them by geometrical factors. 

The finite- difference equations were derived by the application of the 
principle of conservation of mass and momentum to the grid cells. The 
transport processes for each cell are convection and diffusion. Moreover, 
there may be a momentum or mass source within the cell. The mass flux 
across a cell boundary was computed exactly from the scalar product of the 
velocity vector and the vector representing the area of the cell face. 
Note that this can be written out in terms of velocity resolutes and 
geometrical factors including angles between cell faces. In the 
calculation of convection across a cell face, special attention was given 
to the change of the orientation of the coordinate axes from cell to cell 
and the curvature of a cell face. These resulted in extra terras in the 
calculation of convection. However, the representation of convection was 
exact and did not involve any approximation due to the non- orthogonality of 
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the grids. 


The calculation of diffusion is somewhat more complicated than 
convection. The diffusion flux was calculated assuming the coordinate 
system to be locally orthogonal. This obviously neglects cell curvature 
and non- orthogonal orientation and may incorporate a substantial amount of 
error where the process is primarily diffusive. However, in the thin film 
calculation this approximation should not introduce severe inaccuracies, 
particularly when the film enters the control volume with a reasonably high 
velocity. 

The relative importance of convection and diffusion at each cell was 
determined from the magnitude of the local Peclet number. A hybrid 
difference scheme demonstrated by Patankar (1980) was used. The 
calculation of the momentum source due to the pressure gradient and that 
due to the gravitational body force could be accomplished without any 
approximation for non- orthogonality. 

The grid generation was achieved in two steps. First, the grid cells 
were formed by algebraic interpolation between the boundary points. This 
provided an approximately equal volume for each control cell. The 
boundaries for the interior cells were then smoothed to make the cell faces 
more orthogonal to each other. This was achieved by solving the Laplace 
equation for grid geometry. This latter operation resulted in a better 
representation of diffusion in the flow field and more accurate 
computations. The details of the formulation in a body- fitted coordinate 
system and the generation of grid cells are described in the work by 
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Spalding et al. (1986). 


The number of grids in the z-y plane was determined from a series of 
computations with different grid sizes. For radial flow, it was found that 
the free surface profile along with other computed quantities do not change 
any further if the computational mesh is refined beyond 50 x 25 grids in 
the z-y direction. Computations with 50 x 25 and those with 55 x 27 grids 
yielded identical results. For plane flow 40 x 20 grids in the z-y plane 
was found to be adequate, which precisely predicted the friction 
coefficient and velocity profile. in a fully developed falling film flow. 
Therefore, all computations were carried out using 40 x 20 grids for plane 
flow and 50 x 25 grids for radial flow. 

The flow field was solved by using the SIMPLEST algorithm as discussed 
by Spalding (1980). One special feature of this algorithm is that in the 
discretized form of the momentum equation, the convection terms are lumped 
together with the source term. This results in a faster convergence for 
some flow conditions. The algorithm works in an iterative manner where the 
continuity equation is transformed and used as a pressure correction 
equation. The computation starts by guessing a pressure field. This is 
used to determine velocity components from their corresponding momentum 
equations. The modified continuity equation is then used to determine the 
amount of pressure correction. The guessed pressure, the amount of 
pressure correction and the solution corresponding to the momentum 
equations are then assembled together to give the flow rate and pressure 
field for that step. The new pressure serves as a guess for the next step. 
The solution proceeds until the normalized residual for each equation was 
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approximately 10’°. The above finite- difference formulation and solution 
» 

procedure is incorporated in the computer program PHOENICS that was used in 
the present study. 

In the free surface flow discussed here, both the zero- shear condition 
and the p =0 condition at the free surface need to be satisfied. These 
two conditions cannot be simultaneously given at a boundary with the 
existing program. On the other hand, the free surface geometry, which is 
unknown in the problem has to be given before solving the flow field by a 
finite- difference method. To avqid this difficulty, an iteration scheme 
has been adapted as described below. 

(1) Guess a free surface height distribution. One may use the 
one- dimensional solution for a good start. 

(2) Solve the flow field completely for that distribution using the 
zero- shear condition on the free surface boundary. 

(c) Find the pressure distribution on the free surface and calculate 
its deviation from an ideal zero- pressure free surface. The 
measure used here is the normalized RMS (root- mean- square) error 
referenced to the initial total head. It is defined as 

Normalized RMS Error = 

where p^ is the free surface pressure at the kth node, and n is 
the number of nodes adjacent to the free surface. 
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(d) Calculate and reduce the RMS error on the free surface by 
successive alteration of the surface height distribution. 

(e) The results with the minimum error give the required final 
solution. 

The optimization technique used here is known as the exhaustive search 
method. In this method, a general form of the equation representing the 
free surface is assumed. The equation contains a number of arbitrary 
parameters depending on the desired degrees of freedom. To start the 
process, the parameters are given values either from experience (i.e., the 
one- dimensional solution) or just guessed. Only one parameter is then 
changed while the other are held constant and the trend of the RMS error is 
observed. That parameter is changed continuously in successive steps until 
a minimum RMS error is obtained. The parameter is then kept at its optimum 
value while the other parameters are changed one by one following the same 
procedure. Once the first round of optimization is complete, the first 
parameter is changed again in either direction to see whether the error 
increases or decreases. The process is continued with the other parameters 
until an absolute minimum for a combination of parameters is obtained. The 
accuracy of this process depends somewhat on the assumption of the form of 
free surface. In the results presented here, hyperbolas with two or more 
degrees of freedom (Table 3.2) were used to represent the computational 
domain. 

The numerical algorithm had to be modified slightly when a jump was 
present in the flow field. This happened in the case of horizontal film 
flow in the presence of gravity. Dnder some flow conditions, both 
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supercritical and subcritical flows are present, separated by a hydraulic 

i 

jump. When a jump was present, the two regimes of flow ( i.e., 
supercritical and subcritical) had to be computed separately and the 
conditions were matched at the jump interface. It involved the following 
operations. 


(a) Compute the supercritical film height by imposing the inlet film 
height and velocity and by optimizing the free surface profile to give 
the minimum RMS error in the pressure. 

(b) Calculate the jump height from, the supercritical film height using the 
equation 


S 2 

77 




1 + 9.6 Fr 


1 


(3.3) 


where subscript * 1’ denotes the supercritical condition before the 
jump and ’2’ denotes the condition following the jump. This 
relationship can be derived from the momentum balance at the jump 
while the velocity profile before and after the jump are assumed to be 
parabolic in nature. 

(c) Compute the subcritical film height by imposing the flow rate and exit 
film height corresponding to Fr = 1 and optimizing the free surface 
profile to give the minimum RMS error for the pressure. 

(d) The supercritical and subcritical flows are calculated for the same 
flow rate, but for supercritical flow the inlet Froude number is 
fixed, whereas for subcritical flow the exit Froude number is fixed. 

(e) The intersection of the jump height and the subcritical film height 
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determines the location of the jump. Th^ film height before the jump 
is given by the supercritical height whereas after the jump the film 
height is given by the subcritical height. 

In order to increase the computational accuracy for both subcritical 
and supercritical flows, the flow involving a hydraulic jump was always 
assumed to have a parabolic velocity profile at the inlet. Moreover, the 
flow was well-established by evolving the flow with a few extra nodes at 
the upstream side of the computational domain. 

3.5 RESULTS AMD DISCUSSION 
3.5.1 Falling Film 

The flow of a plane film along a vertical wall under the influence of 
gravity is a classical fluid mechanics problem where an analytical solution 
is available for fully- developed laminar flow. In the fully- developed 
region, the film height remains constant and the velocity profile has a 
parabolic appearance. Numerical computations using the present methodology 
were performed for a film which has already reached the fully- developed 
condition. A parabolic velocity profile with the same shape as given by 
the analytical solution was used for the incoming fluid. The Reynolds 
number for the film was Re = 12.5. It was found that for the entire 
domain, the velocity profile across the film remains about the same and the 
friction coefficient was equal to that of the analytical solution. In the 
present investigation, the friction coefficient is defined in terms of 
local average velocity of the fluid across the film in contrast to the 
inlet velocity. This definition is applicable for both plane and radial 
flows. For plane flow with constant thickness, V is constant because of 
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continuity. 


The distribution of pressure on the free surface was also computed and 
the normalized RMS error was determined to be 0.011. This non- zero value 
for the error may be attributed to the inaccuracy associated with the 
representation of the flow field by finite- difference equations. This 
error can be further reduced by more grid points and more degrees of 
freedom in the free surface equation. 

The developing flow of a falling film when introduced at a height 
equal to, above, or below the equilibrium height was also investigated. 
The flow conditions are summarized in Table 3.2 (cases IP) and the results 
are shown in Table 3.3 and Figs. 3.3 and 3.4. Vhen the film enters the 
control volume at the equilibrium height, the height remains the same and 
the development of the velocity profile from uniform to parabolic occurs as 
the flow moves downstream. Vhen the film enters with a height other than 
the equilibrium height, a gradual adjustment of the height takes place 
until the flow reaches the equilibrium height. The adjustment of the free 
surface and the development of the velocity profile occur simultaneously in 
this flow. To model the free surface, a height distribution of the form 
given in Table 3.2 (cases 1PB and 1PC) was assumed. The first part of this 
distribution (£<£.) provides the variation of the free surface height in 
the developing flow region and the second part (^>^) g ives the height 
after the adjustment is completed. The downstream location where the free 
surface adjustment is complete is <f^. 

For the computation with <5^ n = <5.^, the normalized RMS error on the 
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Table 3.3: Summary of Computed Optimum Profiles and 

Normalized RMS Error i-n Pressure 



Problem 

Optimum Profile 

Normalized 
RMS Error in 
Pressure 

Case 1PA 

(Parabolic Inlet) 

uniform 

0.011 

Case 1PA (Uniform Inlet) 

uniform 

0.073 

Case 1PB 


A = 7.5 

an = 0.2 

= 0.2 

0.024 

Case 1PC 


A '= 3000 
an = 0 . 0445 

^ = 0.05 

0.026 

Case 2P 


A = 10 

an = 0.93 

0.054 

Case 2RA 


A = 1.43 

an = -0.64 

0.039 

Case 2RB 


A = 27.3 

an = 0.72 

0.091 

Case 3PA 


A = 1.47 

an = 1.52 

C = 8.42 

D = 1.04 

B = 0.29 

bn = 1.5 

0.014 

Case 3PB 


A = 1.47 
an = 1.83 
C = 11.11 
D = 1.04 
B = 0.29 
bn = 1.5 

0.025 


Continued on next page 
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Case 3RA 


0.033 


A = 

2.05 

an = 

2.2 

C = 

8.51 

. D = 

1.04 

B = 

0.37 

bn = 

1.4 


Case 3RB A 

an 

C 

D 

B 

bn 


=2.05 0.036 

= 2.7 

= 11.01 
= 1.04 

= 0.37 
= 1.4 
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DIMENSIONLESS VERTICAL DISTANCE, £ 

Fig. 3.3 Film height in a developing falling film flow 



free surface condition is 0.073. The uniform entrance condition provides a 

i 

higher RMS error than the parabolic entrance since a developing flow 
regionis present where the pressure has to conform with the flow 
development. An error of this order may be acceptable for a free surface 
computation since no variation of height is expected to take place when a 
film is introduced at the equilibrium height. In the present 
investigation, the normalized RIS error corresponding to the optimum free 
surface profile did not exceed a limit of 0.1. The distance required for 
flow development was found to be about five times the film thickness in 
this case. A flow is defined. to be fully- developed when the friction 
coefficient is within 2 percent of the final equilibrium value. This 
definition is similar to that given by Kays and Crawford (1980) for 
developing flow in closed conduits. 

As seen in Fig. 3.3, the adjustment of the film height takes places 
for a length of approximately 0.2 L, which in terms of equilibrium height 
came out to be 10 when a film is introduced at a height 20 percent more 
than the equilibrium height. In this situation, the velocity profile also 
becomes fully- developed at the same location. A film introduced at a 
height 20 percent lower than the equilibrium is found to require a shorter 
distance for the adjustment of the free surface and the velocity profile. 

The variation of the friction coefficient along the length of the film 
is shown in Fig. 3.4. The distributions for $ in = and <$ in = 0.8^ are 
very close to each other, but for = 1.2<$^ d the variation is 
significantly different. In the last case, the friction coefficient first 
decreases to a value lower than that for the equilibrium condition and then 
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rises up to the equilibrium value. This is v due to two counter- acting 
phenomena that affect the film in the developing flow region in this case. 
The wall friction propagates outward as the boundary layer develops 
beginning from the entrance point which tends to reduce the velocity of the 
fluid. The thickness of the film, however, decreases and tends to increase 
the fluid velocity due to the area available for the flow. Since the first 
effect starts from the wall, it is more dominant in the earlier part of the 
flow development and then the second effect takes over in the region 
downstream. It can also be mentioned that when the film is introduced 
below the equilibrium height, the increase in the film height and the 
propagation of the shear stress tend to reduce the film velocity, so the 
behavior is not analogous to the case when the film is introduced above the 
equilibrium height. In all situations, a plane falling film eventually 
attains a fully- developed flow. This was confirmed by comparison of the 
velocity profile in a location near the exit. The variation of the 
velocity and the friction coefficient were found to be identical in all 
situations . 

3.5.2 Film Flow Under Zero Gravity 

In the absence of gravity, the orientation of the plate becomes 
immaterial and an identical flow condition is achieved if the plate is 
horizontal, vertical, or inclined. A situation is considered herein where 
the film is introduced at a height equal to the equilibrium fully- developed 
flow in a falling film system as discussed in the previous section. In the 
absence of gravity, the flow is acted on only by viscous and inertial 
forces and the film height is expected to increase downstream. To model 
the free surface, a profile of the form given in Table 3.2 (case 2P) is 
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assumed. The flow conditions are also listed in the same table. The 
» 

results are shown in Table 3.3 and Figs. (3. 5-3. 7). 

Figure 3.5 shows the variation of the dimensionless film height with 
distance, which increases monotonically. The figure also shows the 
analytical solution derived in a previous section. The analytical solution 
requires the specification of a friction coefficient. In the present 
investigation, the friction coefficient was taken from the numerical 
solution instead of assuming it to be constant throughout the region. The 
comparison between the analytical and numerical solutions appears to be 
good in most regions of the flow. 

The variation of the w- component of the velocity at three different 
locations is shown in Fig. 3.6. In contrast to the falling film, the 
velocity changes as the flow moves downstream and does not attain a 
fully- developed situation. The profile, however, becomes approximately 
parabolic in nature downstream from the entrance. A test of the velocity 
profile with a true parabola showed that the maximum deviation is less than 
10 percent. The shear stress exerted by the solid wall and the 
corresponding friction coefficient are plotted in Fig. 3.7. It can be 
noticed that the shear stress decreases continuously as the flow moves 
downstream, whereas the friction coefficient has a minimum at an 
intermediate location and then increases. The largest variations of the 
shear stress occur close to the entrance due to the development of the 
velocity from a uniform to a parabolic profile. After the velocity profile 
is fully developed in shape, the slight reduction in the shear stress is 
due to the deceleration of the flow, which is small compared to the 
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DIMENSIONLESS DISTANCE FROM INLET, £ 

Fig. 3.5 Film height for plane flow under zero gravity 





DIMENSIONLESS DISTANCE FROM INLET, f 

Fig. 3.7 Wall shear stress (o) and Friction coefficient (+) for 

plane flow under zero gravity 


reduction of the average velocity in regions, away from the entrance. This 
results in an increase in the friction coefficient. 

The system where a fluid is introduced at the center of a circular 
horizontal plate and spreads uniformly in all radial directions was also 
studied. Two different inlet Reynolds numbers were chosen. The flow 
parameters corresponding to these cases are shown in Table 3.2 (cases 2RA 
and 2RB) . The surface profile and corresponding normalized RMS error for 
the surface pressure are listed in Table 3.3. 

Figure 3.8 shows the variation of the film height with radial 
distance. Case 2RA corresponds to a higher Reynolds number than that of 
Case 2RB. It appears that the inlet Reynolds number is a very strong 
parameter in determining the behavior of the film as it spreads radially 
under zero gravity. At a high Reynolds number, the film decreases in 
height monotonically because the inertial forces are greater than the 
frictional resistance exerted by the solid wall. When the Reynolds number 
is small, the film enters the control volume with a smaller amount of 
inertia and is easily overpowered by the frictional resistance. The strong 
resistance to the flow causes a rather rapid increase in the film height as 
it flows downstream. The figure also shows the film height distribution 
predicted by the one- dimensional analysis for case 2RA. The analytical 
prediction for this case is found to be comparable with the two-dimensional 
numerical solution. 

The variation of the velocity across the film thickness is shown in 
Fig. 3.9 for two different locations. In both cases, the profile is almost 


116 



DIMENSIONLESS RADIAL DISTANCE, ( 

Fig. 3.8 Film height for radial flow under zero gravity 
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Fig. 3.9 Velocity profile for radial flow under zero gravity 



parabolic (deviation within 2 percent) but its . magnitude decreases 
* 

downstream due to the increase in area available for the flow. The 
friction coefficients for these two cases are plotted in Fig. 3.10. In 
case 2RA, the friction coefficient is smaller in magnitude than that of 
case 2RB because the friction coefficient is a function of the wall shear 
stress and the local velocity. Case 2RA, which corresponds to a higher 
fluid velocity, is expected to encounter more resistance from the wall. 
This is precisely the situation as seen in the computed values of the wall 
shear stress. However, the magnitude of the velocity is a stronger 
parameter in the friction coefficient and causes the coefficient to be 
lower in case 2RA. In the figure it can be also noticed that in case 2RB, 
which corresponds to a smaller flow rate, the friction coefficient starts 
increasing soon after the entrance, attains a maximum value at an 
intermediate location and then decreases further downstream. In contrast 
to case 2RA the velocity profile at the inlet was assumed to be parabolic 
in nature. Therefore, a rapid decrease of friction coefficient near the 
entrance, which is characteristic of a developing flow, is not found here. 
The decrease of the friction coefficient downstream indicates that the 
shear stress also decreases very rapidly in that region. This suggests the 
possibility of a zero- shear condition further downstream. 

3.5.3 Horizontal Film Flow in the Presence of Gravity 

The flow of a plane film and a radially spreading film were also 
investigated for a horizontal orientation of the plate where the gravity 
acts across the thickness of the film instead of in the direction of the 
main flow. The flow conditions used here are listed in Table 3.2 (cases 3) 
and the corresponding free surface geometry is described in Table 3.3. 
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Figures 3.11 and 3.12 present the computational results. 

In the case of horizontal flow, two different flow regimes are 
encountered. The flow is subcritical or supercritical depending on the 
Froude number. Both subcritical and supercritical flows move toward a 
critical condition. Vhen the flow enters the control volume as 
subcritical, it remains subcritical in the entire domain since a transition 
from subcritical to supercritical is not possible. However, if the flow 
enters the control volume as supercritical, it may remain supercritical or 
transform into a subcritical flow depending on the amount of inertial, 
gravitational, and viscous forces it encounters during the flow. Moreover, 
if a transition is present it must happen as a jump where the height before 
and after the jump can be related by the jump condition (eqn. 3.3). In the 
present investigation, the supercritical and subcritical branches of the 
flow were computed separately and the location of the jump was determined 

"tlig jump condition. In the computation of supercritical flow, the 
inlet height and Froude number were appropriately prescribed, whereas for 
subcritical flow the outlet condition was prescribed. 

In Fig. 3.11, it can be noticed that the jump moves downstream with an 
increase in the inlet Froude number. Vhen the inlet Froude number is high, 
the flow can maintain its supercritical status for a longer distance. 
Figure 3.12 shows the subcritical and supercritical solutions and the 
corresponding jump height distribution for the two plane flow cases 
considered here. Vhen the inlet Froude number is increased, the inlet 
height decreases if the flow rate is held constant. This causes the 
supercritical height of the film to decrease and the jump height to 
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DIMENSIONLESS HORIZONTAL DISTANCE, £ 

Rg. 3.12 Supercritical and subcritical film height for plane horizontal flow 



increase. Vhen the Froude number is gradually increased at the inlet, a 
condition is reached when the jump height curve goes above the subcritical 
curve and no intersection of the two curves is found in the computational 
domain. At that situation, no jump is possible in the region under 
consideration and the flow remains supercritical over the entire plate. 
Similarly, when the Froude number is gradually decreased at the inlet, the 
jump height curve goes below the subcritical curve. In this situation, no 
supercritical flow can be sustained and the jump occurs at the entrance. 

Comparing the results of plane flow and radial flow, it can be seen 
that the propagation of the jump location downstream with an increase in 
the Froude number is greater for plane flow than radial flow. In the case 
of radial flow, the area increases downstream and causes the flow to slow 
down. This results in more rapid reduction of the Froude number of the 
flow. Therefore, in a radial system a jump is likely to be present for a 
larger range of inlet Froude numbers. 

3.6 CONCLUSIONS 

A numerical solution procedure for the computation of plane or radial 
free surface thin film flows in a normal or zero gravity environment has 
been developed. An analytical solution was also derived for a 
one- dimensional approximation of the flow. A reasonable agreement between 
the numerical and analytical solutions was obtained for most flow 
configurations considered here. Three different flow systems were studied, 
(a) plane falling film, (b) plane and radial flow under zero gravity and 
(c) plane and radial flow on a horizontal plate in the presence of gravity . 
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A plane falling film, whether introduced at its equilibrium height or 
at a height above or below the equilibrium, eventually attains the 
fully- developed condition. In the fully developed region, the numerical 
velocity profile, which is parabolic, matched exactly with the analytical 
solution. An estimate on the error bound of the free surface pressure was 
developed from the calculation of the developing flow of the falling film 
and it was concluded that a normalized RIS error of less than 0.1 may be 
acceptable. The length of the developing region was found to be small. 
Both the free surface height and velocity profile appeared to arrive at the 
equilibrium condition within 10 times the equilibrium film thickness. The 
length was found to be relatively larger when the film begins with a height 
above the equilibrium height. 

For a plane flow under zero gravity, it was found that the film 
thickness monotonically increases as the flow moves downstream. The 
velocity profile is parabolic except for regions very close to the 
entrance. The shear stress at the wall decreases as the flow moves 
downstream. The friction coefficient was computed in terms of the local 
average velocity and was found to increase after coming to a minimum at an 
intermediate location on the plate. 

For radial flow under zero gravity, it was found that at a higher 
Reynolds number the film decreases monotonically in thickness as it spreads 
downstream. At a smaller Reynolds number, there is a relatively rapid rise 
in height and the flow develops a parabolic velocity profile as it moves 
downstream. 
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For both plane flow and radial flow on a horizontal plate in the 
presence of gravity, two different flow regimes (i.e., supercritical and 
subcritical) are found to be present. The transition of the flow from 
supercritical to subcritical is accompanied by a hydraulic jump. The 
location of the jump moves downstream with an increase in the inlet Froude 
number. 
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Section IV 


COMPUTATION OP HEAT TRANSFER IN TVO- DUENSIONAL THIN Fill FLOV 
4.1 SUMARY 

The numerically computed flow field and heat transfer coefficient are 
presented for the flow of a thin liquid film in the presence or absence of 
a gravitational body force. The flow systems studied here include (1) a 
film falling down a vertical wall, (2) plane and radial film flow at zero 
gravity and (3) plane and radial film flow along a horizontal plate in the 
presence of gravity. The heating conditions include isothermal and 
uniformly heated surfaces. The transport conditions considered at the free 
surface are an adiabatic condition when there is no heat loss from the free 
surface and an evaporative free surface maintained at its saturation 
temperature. The height of the free surface, flow field and heat transfer 
coefficient were found to be strongly affected by the magnitude and 
direction of the gravitational body force. They were also found to depend 
on the Reynolds number and Froude number of the incoming fluid. The flow 
conditions changed continuously downstream except for the falling film, 
where a fully- developed condition was established at some downstream 
distance. In the case of horizontal flow in the presence of gravity, a 
hydraulic jump was found to be present under some flow conditions. 
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4.2 imODUCTIOW 


The heating and evaporation of a thin liquid film are found in many 
engineering processes. The desalination of sea water, the distillation of 
petroleum products and evaporation in a cooling tower are examples. The 
use of evaporative cooling is also becoming common in computer and space 
technology because of the high heat rejection requirement per unit surface 
area. An understanding of the flow and heat transfer in a thin film for 
different magnitudes and orientations of the gravitational body force is 
also essential for an appropriate design of the absorber unit of a heat 
pump absorption system in earth and space, which was the motivation for the 
present study. 

The flow and heat transfer in a falling thin liquid film along a plane 
vertical wall has been studied by many investigators since the turn of this 
century. For steady fully- developed laminar flow, a theoretical solution 
can be derived from a simple balance of momentum and energy (Bird et al. 
(1960) and Edwards et al. (1979)). The film height remains constant and 
the velocity profile across the film becomes parabolic in the fully 
developed region. The friction and heat transfer coefficients become 
constant in this region. 

The analysis of developing flow when a film is introduced at its 
equilibrium height is also available in the literature. Faghri and Payvar 
(1979) presented numerical results for laminar flow of a thin liquid film 
down a vertical wall. Both uniformly heated and isothermal surface 
conditions were considered. The effects of evaporation and gas absorption 
on the free surface were also considered in addition to the simple case of 
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heating, where the free surface can be considered to be adiabatic. Their 
results compared reasonably well with the available experimental data. 

A film falling under the influence of gravity ceases to be laminar and 
constant in thickness when the flow rate is high (Re > 50). Vaves tend to 
appear on the surface and the flow becomes turbulent as the flow rate is 
increased. A number of theoretical as well as experimental studies have 
been performed to understand the flow in wavy- laminar and turbulent 
regions. Seban and Faghri (1978) reviewed the experimental data and 
correlations to find the amount, of heat transfer augmentation due to 
surface waves by comparing with results for constant film thickness. Even 
though the nature of the waves could not be identified in general, the heat 
transfer enhancement due to waves appeared to be significant. The study 
covered the cases of heating, evaporation and gas absorption and identified 
the nature of the augmentation peculiar to each of these cases. In a later 
study, Faghri and Seban (1981) presented a theoretical analysis of wavy 
flow assuming a sinusoidal form of the wave. A numerical treatment of 
turbulent flow was presented by Seban and Faghri (1976) for three different 
turbulence models. The numerical and asymptotic solutions were compared 
with their own data as well as previous experiments by Chun and Seban 
(1971). In all these numerical heat transfer studies related to falling 
films, the height of the free surface was assumed to be constant at its 
equilibrium value except for the wavy flow where a sinusoidal form of the 
wave was assumed. 

The radial spreading of a liquid film over a horizontal plate is also 
an interesting free surface problem. Vatson (1964) presented results of 
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analytical and experimental studies of the. radial- spread of a liquid jet 
impinging on a horizontal plane for laminar and turbulent flows. By using 
the boundary layer approximations for the governing equations, analytical 
solutions were derived using a similarity transformation along with the 
Pohlhausen integral method. The analysis covered the regions where the 
boundary layer thickness is less than the film height and where the film is 
totally engulfed by the boundary layer. The effects of the gravitational 
pressure gradient was discussed. The possibility of a hydraulic jump in 
such a flow was also anticipated. However, the analysis was applicable 
only to supercritical flow before, the jump. An equation was presented to 
predict the jump height for any given location of the jump. The agreement 
between the experimental data and the analysis was satisfactory. In this 
study, the fluid, plate, and the surroundings were maintained at the same 
temperature and no heat transfer was present. 

Another area of thin film research is the spreading of a film under 
the action of a centrifugal force as seen in a rotating system. Sparrow 
and Gregg (1959) developed an analytical solution for condensation of 
saturated vapor on a rotating surface. The complete Navier- Stokes and 
energy equations were simplified to a set of ordinary differential 
equations by using a similarity transformation and then integrated 
numerically. Their results gave the condensate layer thickness and the 
heat transfer coefficient along with temperature and velocity profiles. 
Later, Butuzov and Rifert (1972) performed experiments to verify the 
solution of Sparrow and Gregg (1959). In a more recent study, Butuzov and 
Rifert (1973) presented experimental as well as theoretical results for the 
reverse problem of film evaporation from a rotating disk. 
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In all of the previous studies concerning thin liquid films, the 
investigators have tried to develop analytical models or have taken 
experimental data. Some of these models are quite approximate in nature 
and do not bring out the finer details of the flow field. A numerical 
finite-difference solution was attempted only for a falling film flow where 
thickness is uniform and known ahead of time. In general, a free surface 
flow is difficult to solve by the finite- difference method since the 
surface geometry changes along the path of the flow, is unknown ahead of 
time and cannot be fitted in a regular Cartesian or cylindrical coordinate 
system. Moreover, none of the studies mentioned above considered the flow 
under a reduced or zero gravity, which is expected to be different from the 
flow under normal gravity. A proper understanding of such flows is 
essential in the design of space- based cooling systems. 

In a very recent study, Rahman et al. (1989a) (described in section 
III of this report) have developed a finite- difference solution method 
applicable for fluid mechanics of thin film flows under zero and normal 
gravity. A body- fitted coordinate system was used where the free surface 
was approximated by a curve and iterated for the best possible solution. 
The present study is a continuation of that work where an analysis of the 
heat transfer to a thin film is given. In addition to numerical heat 
transfer results, the present study includes an approximate theoretical 
analysis using the Pohlhausen integral method. The results illustrate the 
effects of gravity on the transport for both plane and radial flows. 

4.3 MATHEMATICAL MODEL 

The flow configurations considered in the present investigation are 
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shown in Fig. 4.1. They can be broadly classified into three groups 
according to the presence and orientation of the gravitational body force: 

(1) Falling film along a vertical wall 

(2) Film flow under zero gravity 

(3) Film flow along a horizontal plate in the presence of gravity. 

The three possible orientations with respect to gravity are denoted by 1, 
2, and 3. Moreover, P and R denote plane and radial flows. 

The first case is a classical problem where the major driving 
mechanism is the gravitational body force. This problem is used to check 
the accuracy of the present numerical scheme. The effects of introducing 
the film at a height other than the equilibrium height for a given flow 
rate will be investigated for this case. 

For film flow under zero gravity, two problems will be considered: 

(2P) Plane film flow under zero gravity 
(2R) Radial film flow under zero gravity 

In the absence of any gravitational body force, the orientation of the 
plate becomes immaterial. The flow remains the same whether the plate is 
vertical, horizontal or inclined. In this case the flow is driven by 
inertia and viscous forces. In the radially spreading flow, the area 
available for the film increases downstream and acts as an added mechanism 
for the reduction of the film velocity. 

In the third group, two problems will again be considered according to 
the flow geometry under consideration: 
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(3P) Plane film flow along a horizontal platq, in the. presence of gravity 

I 

(3R) Radial film flow along a horizontal plate in the presence of gravity 

In the case of horizontal thin film flow where the gravitational body 
force acts across the thickness of the film, different flow regimes (i.e., 
supercritical, subcritical or both) may be present according to the local 
film velocity and height. The transition of the flow from supercritical to 
subcritical takes place through a hydraulic jump. Since a sudden 
transition of the flow takes place across a jump, special analytical and 
computational tools are required to, calculate the flow around this singular 
point. 

Two different heating conditions are considered. They are an 

isothermally heated plate (denoted by symbol T) and a uniform heat flux on 
the plate (denoted by H). The thermal condition on the free surface may be 
different depending on whether or not there is any evaporation on that 
surface. In the case of simple heating with no evaporation, the free 
surface may be assumed to be adiabatic in nature. In the case of 

evaporation, the free surface temperature will be the same as the 

equilibrium temperature corresponding to the ambient vapor pressure. In 
the present study, the evaporation is assumed to be small so that the loss 
of fluid at the free surface is negligible compared to the mainstream flow. 

The equations governing the conservation of mass, momentum and energy 
in a thin film may be written as 
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V • V = 0 


(4.1) 


P = - Vp + p V 2 V + pg (4.2) 

^ = a V 2 T (4.3) 

The fluid is assumed to be Newtonian with constant properties over the 
range of temperatures encountered in the problem. 

The body- fitted coordinate system is shown in Fig. 4.2. The local 

coordinate axes are directed along the lines joining the centers of 
adjacent grid cells. The z-axis is directed in the streamline direction 
and the y-axis across the thickness of the film. The velocity resolutes in 
the z- and y- directions are w and v, respectively. The boundary 
conditions are given by: 


at y = 0: 


v = w = 0 

• T = T , for isothermal wall 
-K = q w , for constant flux wall 


(4.4) 


at y = 6: 


r - 0, v = 0, p=0 

f q = 0, for heating 
^ = ^sat’ ^ or eva P ora ti° n 


(4.5) 


at z = 0 or r 


in 


v = V. for uniform entrance 

(4 6 ) 

C 2 (^) ' (J) 2 ] , for parabolic 

entrance 


w = 1.5 V . 

m 
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? ig. 4.2. The coordinate system on a grid cell 
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(4.T) 


at 2 = 1 or r out : H = 0 

, p = 0, for cases 1 and 2 
I P = PS (6 ' y) > for case 3 

The difference between the static pressure in the flow and the ambient 
pressure is denoted by p. On the free surface this quantity has to be zero 
when the surface tension is negligible. 

4.4 NUMERICAL SOLUTION PROCEDURE 

The governing conservation equations (4.1 - 4.3) along with the 

appropriate boundary conditions (4.4 - 4.7) were solved numerically using a 
finite- difference scheme. Since the free surface geometry cannot be 
handled very well with a regular rectangular or cylindrical coordinate 
system, a boundary- f itted curvilinear coordinate system had to be used. In 
this system, the free surface of the film was used as one of the boundaries 
of the control volume. 

As shown in Fig. 4.2, the local coordinates are defined along lines 
joining adjacent cell centers. The z-axis was taken in the streamline 
direction and y-axis in the direction across the film. The resolutes of 
the velocity vectors in the y- and z- directions are v and w, respectively. 
In general, the coordinate system was non- orthogonal in nature, i.e., the 
faces of a grid cell are not orthogonal to each other. The velocity and 
force vectors were resolved into covariant components and the problem was 
solved in its physical domain. 

The finite- difference equations were derived by the application of the 
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principle of the conservation of mass, momentum, and energy to the grid 

i 

cells. The transport processes for each cell are convection and diffusion. 
The mass flux across a cell boundary was computed exactly from the scalar 
product of the velocity vector and the vector representing the area of the 
cell face. Note that this can be written out in terms of velocity 
resolutes and geometrical factors including angles between cell faces. In 
the calculation of convection across a cell face, special attention was 
given to the change of the orientation of the coordinate axes from cell to 
cell and the curvature of a cell face. This resulted in extra terms in the 

i 

calculation of convection. However, the representation of convection was 
exact and did not involve any approximation due to the non- orthogonality of 
the grids. 

The calculation of diffusion was somewhat more complicated than 
convection. The diffusion flux was calculated assuming the coordinate 
system to be locally orthogonal. This neglected the effects of cell 
curvature and non- orthogonal orientation. However, the effects were small 
in the thin film calculation since the flow is convection dominated, 
particularly when the film entered the control volume with a reasonably 
high velocity. 

The relative importance of convection and diffusion at each cell was 
determined from the magnitude of the local Peclet number. A hybrid 
difference scheme demonstrated by Patankar (1980) was used. The 
calculation of the momentum source due to the pressure gradient and that 
due to the gravitational body force could be accomplished without any 
approximation for non- orthogonality . 
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The grid generation was achieved in two steps. -First, the grid cells 
were formed by algebraic interpolation between the boundary points. This 
provided an approximately equal volume for each control cell. The 
boundaries for the interior cells were then smoothed to make the cell faces 
more orthogonal to each other. This operation resulted in a better 
representation of diffusion in the flow field and more accurate 
computations. A more elaborate discussion of the numerical formulation is 
presented in Rahman et al. (1989a) (or section III of this report) and the 
methodology used here is similar to the work by Galea and Markatos (1987) 
where a body- fitted coordinate system was used to predict fire development 
in an aircraft. 

The number of grids in the z-y plane was determined from a series of 
computations with different grid sizes. For radial flow, it was found that 
the free surface height along with other computed quantities do not change 
any further if the computational mesh is refined beyond 50 x 25 grids in 
z-y direction. For plane flow 40 x 20 grids in the z-y plane was found to 
be adequate, which precisely predicted the friction and heat transfer 
coefficients and the velocity profile in a falling film system. Therefore, 
all of the computations were carried out using 40 x 20 grids for plane flow 
and 50 x 25 grids for radial flow. 

The flow field was solved by using the SIMPLEST algorithm as discussed 
by Spalding (1980). One special feature of this algorithm is that in the 
discretized form of the momentum equation, the convection terms are lumped 
together with the source term. This results in a faster convergence for 
some flow conditions. The algorithm works in an iterative manner where the 
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continuity equation is transformed and nsed as a pressure correction 
equation. The computation starts by guessing a pressure field. This is 
used to determine the velocity components from their corresponding momentum 
equations. The modified continuity equation is then used to determine the 
amount of pressure correction. The guessed pressure, the amount of 
pressure correction and the solution corresponding to the momentum 
equations are then assembled together to give the flow rate and pressure 
field for that step. The new pressure serves as a guess for the next step. 
In each step, once the velocity components are known, the temperature field 
is determined by solving the energy equation (4.3). Since temperature is a 
scalar quantity and its equation is linear, the computation of temperature 
is less involved than the velocity components. The solution proceeds until 
the normalized residual for each equation was approximately 10' 6 . In the 
free surface flow discussed here, both the zero- shear condition and the p = 

0 condition at the free surface need to be satisfied. These two conditions 
cannot be simultaneously given at a boundary with the existing program. On 
the other hand, the free surface geometry, which is unknown in the problem 
has to be given before solving the flow field by a finite- difference 

method. To avoid this difficulty, an iteration scheme has been adapted as 
described below. 

(a) Guess a free surface height distribution. 

(b) Solve the flow and temperature fields completely for that height 
using the zero- shear condition on the free surface boundary. 

(c) Find the pressure distribution on the free surface and calculate 
its deviation from an ideal zero- pressure free surface. The 
measure used here is the normalized RMS (root- mean- square) error 
referenced to the initial total head, which is defined as 
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Normalized BIS Error = 


k=l 


pV - 2 

( > 6 Vt" > 

where is the free surface pressure at the kth node and n is 
the number of nodes adjacent to the free surface. 

(d) Calculate and reduce the RMS error on the free surface by 
successive alteration of the surface height distribution. 

(e) The results with the minimum error give the required final 
solution. 

The optimization technique used here is known as the exhaustive search 
method. In this method, a general form of the equation representing the 
free surface is assumed. The equation containes a number of arbitrary 
parameters depending on the desired degrees of freedom. To start the 
process, the parameters are given values either from experience (i.e., 
analytical solution) or just guessed. Only one parameter is changed 
keeping others constant and the trend of RMS error is observed. That 
parameter is changed continuously in successive steps until a minimum in 
RMS error is obtained. The parameter is then kept at its optimum value and 
the other parameters are changed one by one following the same procedure. 
Once the first round of optimization is complete, the first parameter is 
changed again in either direction to see whether the error increases or 
decreases. The process is continued with the other parameters until an 
absolute minimum for a combination of parameters is obtained. The accuracy 
of this process depends somewhat on the assumption of the form of free 
surface. In the results presented here, hyperbolas with two or more 


degrees of freedom were used to represent the v computational domain. 

The algorithm had to be modified slightly when a jump was present in 
the flow field, as in the case of horizontal film flow in the presence of 
gravity. Under some flow conditions, both supercritical and subcritical 
flows are present, separated by a hydraulic jump. When a jump was present, 
the two regimes of the flow ( i.e., supercritical and subcritical) had to 
be computed separately and the conditions were matched at the jump 
interface. It involved the following operations. 

(a) Compute the supercritical film height by imposing the inlet film 
height and velocity and optimizing the free surface height 
distribution to give the minimum RMS error in the pressure. 

(b) Calculate the jump height from the supercritical film height using the 
equation 



1 

5 



1 + 9.6 Fr 


2 

1 



(4.8) 


where subscript ’1* denotes the supercritical condition before the 
jump and ’2’ denotes the condition following the jump. This 
relationship can be derived from the momentum balance at the jump. 

(c) Compute the subcritical film height by imposing the flow rate and exit 
film height and optimizing the free surface profile to give the 
minimum RMS error in the pressure. 

(d) The supercritical and subcritical flows are calculated for the same 
flow rate, but for supercritical flow the inlet Froude number is 
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fixed, whereas for subcritical flow the exit Froude number is fixed. 

(e) The intersection of the jump height and the subcritical film height 
determines the location of the jump. The film height before the jump 
is given by the supercritical height whereas after the jump the film 
height is given by the subcritical height. 

(f) The heat transfer coefficient for the supercritical part of the flow 
is determined from the supercritical temperature field, where the 
inlet temperature to the control volume corresponds to the actual 
inlet temperature. 

(g) The inlet temperature for the subcritical flow is determined from the 
exit bulk temperature corresponding to the supercritical solution. By 
imposing this condition, the energy balance at the jump is 
automatically attained. The heat transfer coefficient in the 
subcritical flow is then determined from the subcritical temperature 
field. 

In order to increase the computational accuracy for the subcritical 
and supercritical flows, the flow involving a hydraulic jump was assumed to 
have a parabolic velocity profile at the inlet. 

4.5 RESULTS AND DISCUSSION 

4.5.1 Falling Film 

The flow and transport in a falling film is a widely studied problem 
in classical thermo-fluid mechanics. In ful ly- developed laminar flow, an 
analytical solution can be derived for the velocity and temperature fields. 
A numerical computation using the present methodology was performed for a 
film which has already reached the fully- developed condition at the 
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entrance. A parabolic velocity profile with the same shape as given by the 
analytical solution was used for the incoming fluid. The Reynolds number 
for the film was Re = 12.5. It was found that for the entire domain, the 
velocity profile across the film remains about the same and the friction 
coefficient was equal to that of the analytical solution. This shows that 
the physical characteristics of a fully- developed flow is retained in the 
numerical solution. The distribution of pressure on the free surface was 
also computed and the normalized RMS error was determined to be 0.011. 
This non- zero value of the error may be attributed to the inaccuracy 
associated with the representation of the flow field by finite-difference 
equations. This error can be further reduced by more grid points and more 
degrees of freedom in the free surface equation. 

The developing flow of a falling film when introduced at a height 
equal to, above, or below the equilibrium height was also investigated. 
The flow conditions are summarized in Table 4.1 and the results are shown 
in Figs. (4.3 - 4.5). When the film enters the control volume at the 
equilibrium height, the height remains the same and the development of the 
velocity profile from uniform to parabolic occurs as the flow moves 
downstream. When the film enters with a height other than the equilibrium 
height, a gradual adjustment of height takes place until the flow reaches 
the equilibrium height. The adjustment of the free surface and the 
development of the velocity and temperature profiles occur simultaneously 
in this flow. 

For the computation where S- ~ S {d , the normalized RMS error on the 
free surface pressure was found to be 0.073. The uniform entrance 
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Table 4.1: Flow Parameters and Free Surface Geometry 


General data: 


Case 1PA: 


Case 1PB: 


Case 1PC: 


Case 2P: 


T. = T , = 10°C 
m sat 

T = 30°C, for isothermal wall 
w ’ 

2 

q = lOOOV/m , for constant flux wall 


S. = Sr. (= 0.000595 m) 
in ia v 

Re =50, Pr = 7, L = 0.03 m 
Free surface: A = A^ 



1,2 ^fd 


Re = 50, Pr 


Free surface: 


7, L = 0.03 m 

1 s (1 + 7.5 O' 0 ’ 2 ’ f ° r $ < °' 2 
A = A fd , for ( > 0.2 


6- = 0.8 S fA 

in fd 

Re =50, Pr 


Free surface 


M 


7, L = 0.03 m 
A = (1 + 3000 £) 0 * 0445 ’ for 
A = A,, , for £ > 0.05 


i < 0.05 


S- = 0.000595 m 
m 

Re = 50, Pr = 7, L = 0.03 m 
Free surface: A = (1 + 10 


(Continued on next page) 
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= 0.1953 m 


Case 2R: <^ n = 0.005 m, 

Re = 1616, Pr = 7, r ifl = 0.0508 m, r Qut 
Free surface: A = (1 + 1.43 £) 


Case 3PA: 



5.7, Fr 


out 


1.0 


Re = 355, Pr = 7, L = 0.1445 m 


Free surface: 

, A = (1 + 1.47 O 1 ’ 52 , supercritical 
{ A = 8.42 [2- (1.04 + 0.29 1 * 5 ] , subcritical (( < 1) 


Case 3PB: 


Fr. 

= 8.6, Fr 

m 


Re = 

355, Pr = 

Free 

surface: 


out 

7, 


= 1.0 

L = 0.1445 m 


1 ft*} 

A = (1 + 1.47 £) , supercritical 

A = 11.11 [2 - (1.04 + 0.29 O 1 - 5 ], subcritical (<f < 1) 


Case 3RA: 


Re in = 202, Pr = 7, r in = 0.0639 m, r Qut = 0.1953 m 
Free surface: 

, A = (1 + 2.05 0 2 ' 2 > supercritical 

{ A = 8.51 [2 - (1.04 + 0.37 £) 1,4 ]i subcritical (( < 1) 


Case 3RB: Fr in = 11.0, Fr Qut = 1.0 

Re in = 202, Pr = 7, r in = 0.0639 m, r Qut = 0.1953 in 

Free surface: 

o 7 

, A = (1 + 2.05 () ' , supercritical 

{ A = 11.01 [2 - (1.04 - 0.37 O 1 * 4 ]* subcritical {( < 1) 
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Fig. 4.4 Nusselt number in falling film flow adjacent to an isothermal surface 










Fig. 4.5 Nusselt number in falling film flow adjacent to a uniformly heated surface 


condition has a higher RMS error than the parabolic entrance since a 
developing flow region is present where the pressure has to conform with 
the flow development. An error of this order may be acceptable for a free 
surface computation since no variation of height is expected to take place 
when a film is introduced at the equilibrium height. In addition to 
obtaining the height distribution with the minimum error, the normalized 
RMS error corresponding to the optimum free surface profile does not exceed 
a limit of 0.1. The distance required for flow development was found to be 
about five times the film thickness in this case. A flow is defined to be 
fully- developed when the wall shear -stress is within 2 percent of the final 
equilibrium value. This definition is similar to that given by Kays and 
Crawford (1980) for developing flow in closed conduits. 

The distance required for the development of the temperature profile 
was found to be somewhat larger than the flow field. For an isothermally 
heated wall, the heat transfer coefficient reached within 2 percent of its 
final equilibrium value at about 206 fd for heating and 40£ fd when 
evaporation was present on the free surface. In the case of evaporation, 
the fluid was assumed to enter the control volume at its saturation 
temperature. A flow with a uniformly heated wall required a somewhat 
longer entry length, which were 256 {d and 45£ fd for heating and 
evaporation, respectively. It can be noticed that pure heating has a 
smaller entrance length than that of evaporation. In both cases, the heat 
transfer coefficient is very high near the entrance and reduces downstream. 
Since heating has a higher heat transfer coefficient, it requires a smaller 
length for development than the length required for evaporation. It is 
also worth mentioning that the definition of the heat transfer coefficient 
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is different in the two cases as given by 
« 

_ r q w ( T v ' V’ for heatin § 

h " ( q w (T w - T gat ), for evaporation ( 4<9 ) 

Figure 4.3 shows the variation of the film height along the plate. 
The adjustment of the film height takes places for a length of 
approximately when the film is introduced at a height of 20 percent 

more than the equilibrium height. In this situation the velocity profile 
also becomes fully- developed at the same location. A film introduced at a 
height of 20 percent lower than the equilibrium is found to take a shorter 
distance for the adjustment of the free surface and the velocity profile. 

The variation of the Nusselt number along the length of the film are 
shown in Figs. 4.4 and 4.5 for isothermal and uniformly heated surface 
conditions, respectively. The Nusselt number here has been defined as 

»« = r ( jr ) 1/3 (4-10) 

This definition is more widely used in falling film literature where 

2 i/3 

(u j g) ' has been used as a length scale. Another way of defining the 
Nusselt number is in terms of the film thickness as given by 

«»* = r (4- ii) 

This is more universal and applicable for zero gravity flows which will be 
discussed in a subsequent section. 
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In all cases, the Nusselt number is very^ high near the entrance. This 

is because the thermal boundary layer starts developing from the entrance 

point on the plate. It decreases downstream and eventually attains an 

asymptotic value when the flow is fully developed. The distribution for 

5^ n = 5^ and = 0.8 5^ coincide with each other and differ by a small 

amount from the distribution for 6 ■ = 1.2 Sr,. A smaller heat transfer 

in id 

coefficient near the entrance when 6- = 1.25,, is due to the fact that 

m Id 

local film thickness there is larger and offers a greater resistance to 

heat transfer. Moreover, the average fluid velocity there is smaller for 

this case. It also appears that the heat transfer coefficient becomes 

uniform downstream whether the film is introduced at a larger or smaller 

height than the equilibrium. This is analogous to the previous observation 

that the film height also attains the equilibrium value after some distance 

* 

from the inlet. The Nusselt number in terms of film thickness, Nu , was 
also computed. It was found that in the case of evaporation, the Nusselt 

number asymptotically approaches unity, the value predicted by the analysis 
of fully developed flow. 

4.5.2 Film Flow Under Zero Gravity 

In a gravity- free environment, the orientation of the plate becomes 
immaterial and an identical flow condition is achieved whether the plate is 
horizontal, vertical, or inclined. A plane flow is first considered where 
the film is introduced at a height equal to the equilibrium fully- developed 
flow in a falling film system as discussed in the previous section. The 
free surface is modeled by assuming a height distribution of the form given 
in Table 4.1. The constants appearing in this equation are optimum values 
determined by the minimization of the error of the free surface pressure. 
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The results are shown in Figs. (4. 6-4. 8). 

Figure 4.6 shows the variation of the film height along the plate, 
which increases monotonically . In the absence of gravity, the flow is 
acted on only by viscous and inertial forces. The effect of inertia is 
maximum at the entrance and decreases downstream. The resistance from the 
wall acts to slow down the fluid, which results in an increment in the 
height in order to preserve the continuity of the flow. The variation of 
the Nusselt number for different boundary conditions for this case are 
shown in Figs. 4.7 and 4.8. The shape of the curves appear to be similar 
for pure heating and evaporation, but differ somewhat between the 
isothermally heated and the uniform heat flux surface conditions. 

* 

The Nusselt number considered here is Nu , which depends on two 

parameters, namely, the film height and the heat transfer coefficient. As 

expected in any developing flow, the heat transfer coefficient is maximum 

at the entrance and decreases downstream. The film height, however, has an 

opposite trend. The net result is the variation presented in Figs. 4.7 and 

4.8. The sudden drop of the Nusselt number near the entrance is due to the 

rapid change of the heat transfer coefficient as the thermal boundary layer 

develops from the leading edge at the entrance section. Except for this 

region, the Nusselt number has a small variation if the plate is 

isothermal. Moreover, when there is evaporation over an isothermal 
* 

surface, Nu is about unity. This shows that heat transfer process becomes 

primarily diffusive after a short distance from entrance. When the plate 

* 

is uniformly heated, the Nu gradually increases downstream after the 
sudden drop close to the entrance. In this case, the reduction of the heat 
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transfer coefficient with distance is relatively smaller and easily 
» 

overpowered by the increment in the film height. Analogous to the falling 
film system, the Nusselt number for heating is always found to be larger 
than that for evaporation. 

The situation where a fluid is introduced at the center of a circular 
plate and spreads uniformly in all radial directions is now considered. In 
a radially spreading flow, the area available for the fluid increases 
downstream and acts as an added mechanism to control the film height and 
the associated transport phenomena.. Since spreading tends to reduce the 
film thickness and friction tends to increase it, the height of the film 
may increase or decrease downstream depending on the flow rate and inlet 
height. The behavior of the flow may be characterized by the inlet 
Reynolds number. The optimum film height distribution for Re in = 404 is 

given in Table 4.1 and is graphically shown in Fig. 4.6. For this 
particular flow condition, the film height decreases monotonically 
downstream. 

The Nusselt number variations for radial flow considered here are 

shown in Figs. 4.9 and 4.10. Analogous to the plane flow case, there is a 

* 

sudden drop of Nu close to the entrance. For a uniformly heated surface, 
the Nusselt number keeps decreasing downstream for the heating and 
evaporation situations. In most developing flows, the heat transfer 
coefficient is maximum at the entrance and decreases with the growth of 
thermal boundary layer. In this case, the film height also decreases 
downstream, so this trend in the Nusselt number is expected. For an 
isothermal surface, however, the Nusselt number drops to a minimum at an 
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intermediate location on the plate and then gradually increases as the flow 
moves further downstream. This trend also indicates that, for an 
isothermal surface, the heat transfer coefficient has to have a minimum at 
an intermediate location on the plate. In the results of the heat transfer 
coefficient, it was found that a minimum is indeed present for both 
isothermal and uniformly heated surfaces. This phenomenon can be explained 
by the fact that in a radially spreading flow with the rapid reduction of 
the film thickness, the thermal boundary layer engulfs the entire film 
downstream from entrance for a moderate Prandtl number fluid considered 
here. Therefore, with a further reduction of film thickness, the thermal 
boundary layer thickness also decreases and results in a reduction of the 
resistance to heat transfer from the wall to the fluid. 

It is also worth mentioning that for both plane and radial flows at 
zero gravity, the flow field never came to an equilibrium fully- developed 
condition. The velocity profile across the film thickness, however, 
appears to be parabolic at most locations of the flow. 

Horizontal F ilm Flow in the Presence of Gravity 
The flow of a plane film and a radially spreading film were also 
investigated for a horizontal orientation of the plate where the gravity 
acts across the thickness of the film instead of in the direction of the 
mam flow. The flow conditions used here and the corresponding free 
surface height distribution are listed in Table 4.1. Figures (4.11-4.16) 
present the computational results. 

Two different flow regimes are encountered in the case of horizontal 
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Fig. 4.11 Film height for plane horizontal flow in the presence 


Fig. 4.12 Nusselt number in plane horizontal flow 

isothermal surface 



Fig. 4.13 Nusselt number in plane horizontal flow on a uniformly heated surface 
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Fig. 4.15 Nusselt number in radial horizontal flow on an 

isothermal surface 
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Fig. 4.16 Nusselt number in radial horizontal flow on a 

uniformly heated surface 


flow. The flow is either subcritical or -supercritical depending on the 
Froude number. The subcritical and supercritical flows move towards a 
critical condition. Vhen the flow enters the control volume as 

subcritical, it remains subcritical in the entire domain since a transition 
from subcritical to supercritical is not possible. However, if the flow 
enters the control volume as supercritical, it may remain supercritical or 
transform into a subcritical flow depending on the amount of inertial, 
gravitational, and viscous forces it encounters during the flow. Moreover, 
if a transition is present it must happen as a jump where the height before 
and after the jump can be related. by the jump condition [equation (4.8)]. 
In the present investigation, the supercritical and subcritical stems of 
the flow were computed separately and the location of the jump was 
determined via the jump condition. In the computation of supercritical 
flow, the inlet height and Froude number were appropriately satisfied, 
whereas for subcritical flow, the outlet condition was satisfied. 

In Figs. 4.11 and 4.14, it can be noticed that the jump moves 
downstream with an increase in the inlet Froude number. Vhen the inlet 
Froude number is high, the flow can maintain its supercritical status for a 
longer distance. It can also be seen that the propagation of the jump 
location downstream with an increase in the Froude number is greater for 
the plane flow than the radial flow. In the case of radial flow, the area 
increases downstream and causes the flow to slow down. This results in 
more rapid reduction of the Froude number of the flow. Therefore, in a 
radial system a jump is likely to be present for a larger range of inlet 
Froude numbers. 
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The distribution of the heat transfer coefficient for different 

thermal boundary conditions are plotted in Figs. 4.12 and 4.13 for plane 

flow and in Figs. 4.15 and 4.16 for radial flow. It appears that the 

Nusselt number decreases rapidly close to the entrance as the thermal 

boundary layer develops. Also, the Nusselt number increases at the 

location of the jump. A sudden transition of the flow happens at the jump 

and the thermal boundary layer starts growing again after the jump. In the 

present study, the heat transfer in the supercritical and subcritical flows 

were computed separately by matching the temperature and heat flux at the 

jump location. Therefore, the .heat transfer coefficient at the jump 

location may not be correct. However, the trend that the heat transfer 

coefficient increases at the jump is consistent with the physical behavior 

of the flow. Except for the jump and leading edge of the boundary layer, 

the Nu decreases slowly with the horizontal distance in the supercritical 

part of the flow and remains approximately constant in the subcritical stem 

of the flow. In Figs. 4.11 and 4.14, it is noticed that for both plane and 

radial flows, in supercritical regime the film height increases with 

distance, whereas in the subcritical regime it decreases with distance. 

Vhen the height of the free surface increases and the fluid velocity 

decreases, the heat transfer from the plate to the fluid is expected to 

decrease downstream. In subcritical flow, however, the fluid velocity 

itself is smaller and a balance is established between the thickness of the 

thermal layer and the velocity of the fluid near the wall. Figure 4.13 

* 

also shows a plot of Nu for the case of heating on a horizontal uniform 
flux surface. It can be noticed that the increment of Nu at the jump is 
much larger than that corresponding to Nu. Since the film height and heat 
transfer coefficient both increase at the jump, this rise in Nu* is quite 
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expected. The figure also shows that the numerical value of Nu approaches 

i 

the value predicted by theoretical solution. 

4.6 CONCLUSIONS 

Numerically computed results for heat transfer in a thin liquid film 
are presented. The flow situations considered are: (a) plane falling film, 
(b) plane and radial flow under zero gravity and (c) plane and radial flow 
on a horizontal plate in the presence of gravity. In all cases, two 
thermal boundary conditions at the wall, namely, isothermal and constant 
heat flux, and two free surface thermal boundary conditions, namely, 
adiabatic and evaporative are considered. A curvilinear body- fitted 
coordinate system is used to handle the irregular flow geometry. 

It was found that in a falling film system, an equilibrium 
fully- developed region is present where the gravitational body force is 
balanced by the viscous shear force at the wall. The film, whether 
introduced at the equilibrium height or at a height above or below the 
equilibrium, eventually attains this fully- developed condition. In the 
fully developed region, the numerical velocity profile matched exactly with 
the analytical solution irrespective of the entrance condition of the film. 
The heat transfer coefficient also attained its fully- developed value and 
matched well with the analytical solution for all of the entrance 
conditions . 

The length required for the development of the flow was found to be 
small. Both the free surface height and the velocity profile appeared to 
arrive at the equilibrium condition within 10 times the equilibrium film 
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thickness. The development length was found to be relatively larger for 
heat transfer. A length of 20-25 times the film thickness was required for 
the development of the equilibrium heat transfer conditions when the film 
was just heated and no evaporation was present on the free surface. In the 
presence of evaporation, the development length was even longer, being 
about 40-45 times the film thickness. There was a small difference in the 
heat transfer coefficient and the development length for different entry 
conditions considered here. 

For a plane flow under zero gravity, it was found that the film 
thickness monotonically increases as the flow moves downstream. For an 
isothermally heated wall, the heat transfer coefficient gradually decreased 
downstream. The Nu* in this situation was approximately constant except 
for regions very close to the entrance. For a uniformly heated wall, the 
Nu reached a minimum at an intermediate location close to the entrance and 
then increased slowly downstream. 

In the case of radial flow under zero gravity, it was found that for 
the flow rate considered here the film thickness decreases monotonically as 
it spreads downstream. In this situation, the Nusselt number was also 
found to decrease monotonically when the plate is uniformly heated. For an 
isothermal wall, however, the Nusselt number reached a minimum and then 
slowly increased downstream. This behavior of the heat transfer is related 
to the thinning of the film and growth of thermal boundary layer in a 
diffusion- dominated heat transfer process. 

For both plane flow and radial flow on a horizontal plate in the 
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presence of gravity, two different flow regimes (i.e., supercritical and 
subcritical) are found to be present. The transition of the flow from 
supercritical to subcritical is accompanied by a hydraulic jump. The 
location of the jump moves downstream with an increase in the inlet Froude 
number. The heat transfer coefficient is discontinuous at the jump. In 
the supercritical region, it decreases gradually downstream, whereas in the 
subcritical region, it remains approximately uniform at most locations on 
the plate. 
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Section V 


DEVELOPMENT OF A NEV COMPUTATIONAL METHODOLOGY FOR THE 
FREE SURFACE FLOWS USING A PERMEABLE WALL 


5.1 SUMMARY 

A new computational procedure for determining the structure of the 
free surface flow of a thin liquid film is presented. The iterative method 
assumes the free surface to be a porous wall where transpiration through 
the wall is allowed while maintaining it at a constant pressure condition 
with no shear stress. The Eulerian computation uses a body- fitting 
coordinate system and an iteration procedure where successive improvements 
of the free surface geometry is attained from the velocity components on 
the free surface. In the final iteration, the transpiration becomes 
negligibly small and thereby the free surface forms a streamline. This new 
algorithm has the advantage over existing computational methods in that a 
complete two-dimensional solution of the flow field and heat transfer 
coefficient can be obtained and can be applied to complex flow problems 
like a hydraulic jump. The computed results include plane and radial flows 
involving a hydraulic jump and those flows at zero gravity where no jump 
can be present. The details of the flow structure, the friction 
coefficient and the heat transfer coefficient are presented. 
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5.2 imODUCTIDff 


Free surface flows are encountered in a vide variety of physical 
phenomena. The flow in a river or water channel, the flow of a thin film 
during evaporation or condensation, and the impingement and spreading of a 
liquid jet on a solid wall are a few examples. All of these phenomena 
share a common feature in that the domain of interest has an unknown 
boundary, the structure of which depends on the flow and the ambient 
conditions. The complexity of the subject is associated with the knowledge 
of the free surface location in conjunction with the satisfaction of the 
appropriate boundary conditions and governing equations. These flows are 
also heavily influenced by the magnitude and direction of the gravitational 
body force. The motivation of this study is to develop a comprehensive 
computational tool for free surface flows both on earth and in space. 

Like any external flow adjacent to a solid surface, the viscous 
effects in a free surface flow are also confined to a boundary layer 
adjacent to the solid boundary. Therefore, inviscid flow models can be 
used for flows involving a depth much larger than the boundary layer 
thickness. The open channel flows in a river fall into this category. An 
extensive review of numerical methods for potential flow involving a free 
surface is presented by Yeung (1982). His review covered both linear and 
non-linear flow problems and three major categories of methods: 
finite- difference, finite element, and boundary integral equations. The 
possibility of using a boundary- fitted coordinate system for the 
representation of free surfaces was also discussed. He concluded that 
hybrid methods based on matching an interior numerical solution with an 
exterior analytical representation is the most rational procedure. 
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The potential flow solution becomes inadequate when viscous effects 
are important for the flow under consideration, and the solution of the 
complete Navier- Stokes equation becomes necessary. The numerical 
techniques presently in use for this class of problems can be divided into 
three categories, namely, (a) surface- tracking methods, (b) volume- 
tracking methods and (c) moving- grid methods. 

In surf ace- tracking methods, the interface is specified by an ordered 
set of imaginary points; between these points its position is approximated 
by an interpolant, which is usually a piecewise polynomial. This 
time- dependent interface divides the flow domain into connected regions in 
which different fluids exist. One of the variants of the surf ace- tracking 
method is the height- function method used by Hirt et al. (1975). The free 
surface was represented by its distance from a fixed surface. This 
distance (height) changed with time and was governed by an equation 
expressing the fact that the interface must move with the local flow field. 
A more commonly used technique to track a surface is to track a string of 
imaginary particles spread along the density interface. Awn (1979) applied 
this method for solving flame propagation. Jun (1986) developed an 
improved version of the interface- tracking method where the location of the 
interface was solved by a scalar equation. He applied his technique to 
compute the sloshing motion of a liquid in a tilted tank, as well as the 
drainage of a pipe under normal and zero gravity conditions. 

In a volume- tracking method, the interface is specified by the common 
boundary of the two regions adjacent to each other. The region is 
identified by its possession of fluid markers of a particular kind. The 


174 



reconstruction of the actual surface shape requires the knowledge of fluid 
markers in the area under consideration. One of the earliest volume- 
tracking methods is the Marker and Cell (MAC) method of Harlow and Velch 
(1965) . Marker particles are scattered initially to identify each fluid 
region in the calculation domain. These particles are transported in a 
Lagrangian manner along with the fluids. Amsden and Harlow (1970) 
developed a simplified MAC method called SMAC, which had some mathematical 
advantages over its predecessor. More recently, Hirt and Nichols (1981) 
designed a different variation of the volume- tracking method called the 
fractional volume of fluid (VOF) method, which requires less storage than 
the MAC method. 

In moving- grid methods, the original grid system can be adjusted to 
approximate the interface. The governing equations are then solved on the 
new distorted grid by treating the grid points on the density interface as 
a moving boundary. In this manner, undesirable numerical mixing between 
the different fluid regions is reduced or avoided. One way to formulate a 
moving grid algorithm is to use a rectangular Cartesian or cylindrical 
coordinate system and at each iteration move the interface to the nearest 
cell wall so that a grid cell always contains one particular fluid. A more 
improved method, as used by Jun (1985), is to partially block a grid cell 
to form the interface within a rectangular Cartesian or cylindrical grid 

cell. 

The numerical methods described above have mostly used a Lagrangian 
approach for solution and dealt primarily with transient problems. 
Moreover, none of these studies have considered low Reynolds number flow 
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like a thin film flow where viscous effects are predominant. In a very 
recent study, Rahman et al. (1989a) (described in section III of this 
report) developed a new moving- grid method, where a boundary- fitted 
coordinate system was used and the free surface was taken as one of the 
boundaries of the computational domain. The primary emphasis of that study 
was to understand low Reynolds number flows as encountered in a thin film, 
both in normal and zero gravity environments. Both plane flow over a 
surface and radially spreading flow were considered. Their results showed 
that the gravitational body force has a strong influence on the structure 
of the free surface. A hydraulic jump was found to be present when a 
gravitational body force acted across the thickness of the film. In that 
situation, two distinct flow regimes were found to be present, which are 
supercritical flow upstream from the jump and subcritical flow downstream 
from the jump. In that study, the free surface was represented by an 
analytical equation with two or more arbitrary constants. The constants 
were optimized with an exhaustive search method where the final solution 
most closely matched the atmospheric (or ambient) pressure distribution on 
the free surface. 

Even though the study by Rahman et al. (1989a) handled most low 
Reynolds number flows fairly well, difficulties were encountered in 
handling flows involving a hydraulic jump, where two regions of the flow 
(supercritical and subcritical) could not be computed as a single- domain 
problem. The computations were done separately for the two regions and the 
solutions were matched at the location of the jump. This procedure, even 
though the location of the jump was correctly predicted, could not account 
for the details of the flow field and heat transfer behavior in the 
vicinity of the jump. An effort to remove this drawback led to the 
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invention of the present algorithm, which can solve the problems involving 
a hydraulic jump as a single domain problem. It is essential to have this 
feature in the calculation of thin film flow over a rotating surface, which 
is one of the motivations for the initiation of the present efforts. It is 
also found that the present algorithm is more universally applicable to any 
free surface problem in the presence or absence of a gravitational body 
force. Unlike any other method presented above, the present procedure 
assumes the free surface to be a permeable wall through which fluid may 
leave or enter the control volume depending on the ambient and local static 
pressures. The location of the surface is improved by successive 
iterations until penetration through the free surface become negligible. 


I 

4 
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5.3 GOVERNING EQUATIONS 

The present numerical methodology will be tested against a number of 
steady free surface problems both in normal and zero gravity environments. 
The problems considered here are schematically shown in Fig. 5.1. They are: 

(a) The two-dimensional flow of a plane thin film over a plate, ' as 
introduced from a pressurized fluid reservoir in a gravity- free 


environment. 

(b) The radially spreading flow of a thin film over a plate as initiated 
by an impinging fluid jet or pressurized reservoir at the center of 
the plate in a gravity- free environment. 

(c) The plane flow over a horizontal plate in the presence of gravity 
where a hydraulic jump is present. 

(d) The radially spreading flow over a horizontal plate in the presence of 
gravity where a hydraulic jump is present in the computation domain. 
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The coordinate system on a grid cell 





Plane flow under zero gravity 



/////////; VV ///V/7 

Radial flow under zero gravity 



Plane hydraulic jump Radial hydraulic jump 


Fig. 5.1. The coordinate and flow systems in the present investigation 
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In all of these problems, the heat transfer behavior in the liquid 
layer will be studied by assuming the solid wall temperature to be higher 
than the temperature at which liquid enters the control volume. The wall 
also remains at an isothermal condition for the entire length of the flow. 
In the case of simple heating with no evaporation, the free surface may be 
considered to be adiabatic in nature. In the case of evaporation, the free 
surface remains at an isothermal condition with a temperature equal to the 
saturation temperature corresponding to the ambient vapor pressure. For 
the problems considered here, the evaporation is assumed to be small so 
that the loss of fluid at the free surface is negligible compared to the 
mainstream flow. 

The equations governing the conservation of mass, momentum and energy 
in a liquid layer involving a Newtonian, constant- property fluid can be 
written as 

V • V = 0 (5.1) ' 

-4 

H = - ivp ♦ </V 2 V ♦ g (5.2) 

Si = OT 2 T (5.3) 

Figure 5.1 also shows the coordinate system used here. The local 
coordinates are directed along the lines joining the centers of the 
adjacent grid cells. The z-axis is directed along the direction of the 
mainstream flow and y-axis is directed across the thickness of the film. 
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The velocity resolutes in the y- and z- directions are v and w, 
respectively. The boundary conditions for the problems considered here are 
given by 


aty=0: v=w=0, T = T y 

(5.4) 



at y = 5: r.n = 0, l«t = 0 


,q = 0, for heating 

(5.5) 

IT = T j., for evaporation 


at z = 0 or r in : w = V in , for uniform entrance 

(5.6) 


v = 1.5 V in [2 (^) - (^) 2 ], for parabolic entrance 
v = 0 


at z = L or r Qut : ^ = 0, P = Pg ($ - y) ( 5 * 7 ) 

Fr = 1, for flows involving a jump 

The normal stress condition on the free surface leads to equations 

balancing the pressure and other stresses including surface tension. 

* 

Scaling these equations, one can show that, for thin film flow where Re 
and Ve are large, the pressure on both sides of the free surface must be 
the same. The quantity p is defined to be the difference between the 
static pressure in the flow and the ambient pressure. In the absence of 
any significant surface tension and tangential stress from the ambient 
fluid, this quantity becomes zero on the free surface. At the outlet end, 
the pressure gradient is hydrostatic in nature. The local parabolic 

approximation in numerical computation eliminate any rigorous specification 
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of v and w at that point. 


5.4 SUmiCAI SQLPTIOW PROCEDURE 

A finite- difference numerical scheme was used for the solution of the 
governing equations (5. 1-5.3) along with the appropriate boundary 
conditions (5.4- 5.7) for the different problems considered here. To 
properly handle the irregular free surface geometry, a curvilinear 
body- fitted coordinate system was used, where the free surface was taken as 
a coordinate surface. 

The computational domain was divided into a non-uniform mesh where 
more grid cells were placed at locations where steeper variations of the 
flow and heat transfer are expected to take place, like the vicinity of a 
hydraulic jump. The cells were formed by an algebraic interpolation 
between boundary points. In general, the cell faces were non- orthogonal in 
nature . 

* 

The local coordinates were defined along lines joining the adjacent 
grid cells as shown in Fig. 5.1. The z-axis is in the streamline direction 
and the y-axis is across the thickness of the fluid layer. The velocity 
vectors were resolved into co- variant components and the problem was solved 
in its physical domain. 

The finite- difference equations were derived by the application of the 
principles of the conservation of mass, momentum and energy to the grid 
cells. The derivations used ideas similar to those presented by Rahman et 
al. (1989a) or section III of this report. The hybrid difference scheme, 
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as demonstrated by Patankar (1980) -as also incorporated in the 
formulation. This basically determined the total contribution (convection 
and diffusion) to a cell from its neighbors as a function of the cell 

Peclet number. 

The flow field .as solved by using the SIIPLEST algorithm as presented 
by Spalding (1980) . One special feature of this algorithm is that in the 
discretized form of the momentum equation, the convection terms are lumped 
together with the source term. This results in a faster convergence than 
its predecessors SIMPLE or SIMPLEM as described in Patankar (1980) . The 
SimEST method also works in a iterative manner .here the continuity 
equation is transformed and used as a pressure correction equation. The 
computation starts by guessing a pressure field. This is used to determine 
the velocity components from their corresponding momentum equations. The 
modified continuity equation is then used to determine the amount of 
pressure correction. The guessed pressure, the amount of pressure 
correction and the solution corresponding to the momentum equations are 
then assembled together to give the flow rate and pressure field for that 
iteration step. The new pressure serves as a guess for the next step. The 
discretized equations for the velocity components .ere solved slab- by- slab 
marching in the z-direction, whereas the pressure and temperature were 
solved by the whole- field method. The convergence of the solution was 
monitored by examining the normalized residuals for each equation and by 
spot- checking values of the pressure, velocity components, and temperature 
at critical locations in the flow. The computation proceeded until the 
magnitude for all of the residuals became negligibly small and the 
spot- checked values attained an invariant condition. The computed results 
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were also tested for grid independence by systematically varying the number 
of grids in both the y- and z- directions. A value of 25 grids was found to 
be adequate in the y- direction, whereas 50 to 150 grids were used in the 
z-direction for the problems considered here. 

In any computational scheme for a free surface flow, the major 
difficulties encountered are in handling the free surface itself. First, 
the geometry of the free surface needs to be specified before solving the 
flow using a finite- difference scheme. On the other hand, the geometry 
itself is dependent on the flow conditions. Moreover, a number of boundary 
conditions needs to be satisfied: (a) the free surface is a streamline (the 
velocity vector on the surface must be tangent to the surface itself at all 
locations) (b) a zero- shear condition must exist on the surface and (c) in 
the absence of any significant surface tension, the static pressure of the 
fluid next to the free surface must be equal to the ambient pressure. All 
of these conditions are satisfied simultaneously for only the correct 
location of the free surface, and therefore some of the conditions ‘are 
violated when an assumed free surface profile is used to initialize the 
solution process. Thus, one has a choice of which boundary condition he 
wants to satisfy exactly. In the present algorithm, the boundary 
conditions (b) and (c) are satisfied exactly, whereas an iteration process 
is introduced to arrive at condition (a). 

It is assumed that the free surface is a permeable wall through which 
fluid particles may leave or enter depending on the pressure on the two 
sides of the surface. The ambient pressure is prescribed. An outflow 
takes place when the static pressure of the fluid inside the control volume 


183 



is higher than the ambient pressure, whereas an inflow occurs when the 
ambient pressure is higher. Since inflow or outflow is permitted, the 
velocity vector at the surface is no longer aligned with it, so an 
iterative scheme can be devised to adjust the surface so that a streamline 
condition on the free surface can be approached. Another criterion to 
follow is to minimize the loss or gain of fluid through the surface. The 
scheme works as follows: 


(1) Prescribe a free surface height distribution. 

(2) Solve the flow field completely for that height using p = 0 as the 

boundary condition on the free surface. 

(3) Find the amount of penetration of the fluid through the surface at all 

locations along the flow. 

(4) Calculate the deviation from the zero- penetration condition. The 
measure used here is the normalized root- mean- square of the 
penetration. 
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(5) Calculate the new free surface height by adding a correction to the 
old height. As shown below, the correction is determined from the 
Cartesian components of the velocity vector on or near the free 
surface by satisfying the condition that the free surface is a 
streamline. The method is analogous to the moving height function 
approach used by Hirt et al. (1975). 
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(6) Continue iterating until the normalized penetration attains a minimum 
value. 

(7) Once the location of the free surface is determined, solve for both 
the flow and temperature fields for that surface profile. 

The correction in the film height is determined from the penetration 
rate at the free surface. The velocity component across the film that is 
associated with the penetration can be estimated to be 

v io SS = T c - as "c ( 5 - 


where v c and w £ are the Cartesian components of the velocity vector at the 
cell next to the free surface. Ve need to drive v^ oss to zero to arrive at 
the streamline condition on the free surface. Therefore, a new height can 
be determined as 




(5.9) 


A central difference is used to calculate the slope of the free surface 
from the height in the previous iteration. This resulted in the equation 



An alternative formulation for the film height correction was also 
derived by using the HacCormack predictor- corrector method, as presented by 
Anderson et al. (1984), which resulted in the formulation: 
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for w c > 0 
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(5.11b) 


The correction scheme needed to be modified slightly at the downstream 
boundary since only the values on one side of the cell were available. The 
upstream boundary did not pose any problem since the film height there 
remained constant. Computations with either equations (5.10) or (5.11) 
were found to be equally good, and lead to the same final solution. 


The procedure required a reasonably good guess of the initial free 
surface height for uniform convergence within a reasonable number of 
iterations. Since penetration is allowed on the free surface, a surface 
profile far from the actual free surface produces outflow on the ordef of 
the mainstream flow and leads to divergence. Therefore, an initial guess 
is limited to a surface which yields penetration smaller in magnitude than 
the mainstream flow. In the calculations presented here, the start-up 
solution was calculated using a ’pressure- optimization method’ as presented 
by Rahman et al. (1989a) or section III of this report. Since that method 
is robust and was available to us, it was used here to obtain faster 
convergence. However, one may start from an approximate analytical 
solution given in section VI or the one- dimensional numerical solution 
presented in section I. However, as expected the convergence time depends 
on the initial condition. Hence, a crude starting solution may take a 
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relatively longer period for convergence. The recommended procedure is to 
use the present method in combination with the ’pressure- optimization 
method’ presented by Rahman et al. (1989a). The former doing the fine 
tuning on the global variations predicted by the latter. 

The correction scheme for the free surface height may require 
under- relaxation to make the algorithm stable and uniformly convergent. 
Also, some smoothing may be required to get around fluctuations that are of 
the same order of the grid size. Vhen the fluctuations were large, a 
polynomial fit of the new height was used. As used by Hirt et al. (1975), 
the smoothing of the free surface from one iteration to the next is 
required to preserve a better continuity of the free surface profile. 

In order to dampen out the effects of flow development, a parabolic 
velocity profile is assumed to be present at the inlet during the iteration 
procedure. Since penetration through the surface is primarily a function 
of the pressure and velocity components, and the energy equation is not 
coupled to the continuity or momentum equations, the flow field is solved 
first as an isothermal problem during the iterations and the energy 
equation is solved only when the free surface location is established. 
This minimizes the computational effort to some extent. 

5.5 DISCUSSION OF RESULTS 

5.5.1 Film Flows in the Absence of Gravity 

Ve now consider the flow of a thin film in a gravity- free environment. 
In the absence of any other body or surface force, the flow is acted on 
only by inertial and viscous forces. Therefore, the orientation of the 
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solid wall, whether vertical, horizontal or inclined, does not play any 
role in the behavior of the flow. Both plane and radial flows are 
considered here. The results are presented in Figs. 5.2- 5.6. 

The predicted variation of the film height for plane film flow under 
zero gravity is presented in Fig. 5.2. The Reynolds number for the flow is 
Re = 50, which is in the laminar flow regime. The properties used are for 
water at 20°C under normal atmospheric conditions. The film height 
increases monotonically with distance. As the fluid moves downstream, it 
decelerates due to friction and consequently the height of the film 
increases to preserve the continuity of the flow. The figure also shows 
the film height predicted in previous investigations. The previous results 
for two-dimensional flow using the pressure optimization method is quite 
close to the present numerical solution. This is expected since both 
studies preserved the two-dimensional behavior of the flow field, and did 
not make any approximations to simplify the problem. One constraint on the 
free surface geometry in the studies of Rahman et al. (1989a) is that .the 
free surface conforms to an algebraic equation with two or three degrees of 
freedom. This constraint is removed in the present study, where the free 
surface can conform to the local flow conditions in a more exact fashion. 
The degree of freedom attained here is the same as the number of grids in 
the main flow direction. Even though an algebraic equation is easier to 
use, the discrete representation of the free surface as done here leads to 
more accurate results. Figure 5.3 shows the distribution of the penetration 
rate through the free surface. Comparing the amount of flow through the 
free surface, it was found that the present free surface has about 307. less 
penetration than the surface predicted by the method in Rahman et al. 
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Fig. 5.2 The film height variation for plane flow under zero gravity 
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Fig. 5.3 Penetration through the free surface for plane flow under zero gravity 
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Computed results for plane film flow under zero gravity 
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(1989a) for the same flow conditions. The average penetration rate through 
the surface was found to be only 0.67. of the entrance flow. This suggests 
that the present methodology is practically applicable, and leads to 
results that are more accurate than the previous two-dimensional method. 

Figure 5.2 also compares the film height with that predicted by Thomas 
et al. (1989) where the variations of the velocity across the thickness of 
the film were integrated out using a parabolic velocity profile and the 
height was explicitly solved by using the MacCormack predictor- corrector 
method. The one- dimensional solution seems to predict a smaller film 
height than either the present or previous two-dimensional numerical 
solution. In addition to neglecting any variation of the velocity across 
the thickness of the film, the one- dimensional solution assumes a friction 
factor based on the boundary layer thickness. These approximations lead to 
results that are somewhat different than the two-dimensional solution 
presented here. Looking at Fig. 5.3, it can be noticed that the surface 
predicted by the method proposed by Thomas et al. (1989) results ' in 
penetrations that are much higher than the solution presented here. 
Comparing the present solution with the analytical results presented in 
section VI, the present film height is found to be lower than that 
predicted by the analytical solution. The analytical solution assumes a 
uniform velocity profile across the film thickness, which is apparently not 
the case in reality. The friction factor used for the final solution, 
however, was the average friction factor obtained from the present 
numerical solution. The above observations suggest that the free surface 
flow of a thin film is two-dimensional in nature and the one- dimensional 
results (analytical or numerical) can be used only as an approximation. 
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The variation of the friction coefficient and the Nusselt number for 
plane flow under zero gravity are presented in Fig. 5.4. The friction 
coefficient decreases rapidly close to the entrance and gradually increases 
thereafter all the way to the exit. The steep drop near the entrance is 
associated with the development of the velocity profile from a uniform 
entrance condition to a parabolic configuration downstream. By running the 
problem with a parabolic entrance flow, it was found that the rapid drop 
near the entrance mostly disappeared, but the friction coefficient after a 
short distance from the entrance coincides with the plot presented here. 
One may also notice that the friction coefficient is calculated in terms of 
the average local velocity, which decreases gradually downstream. 
Observing the numerically calculated values of the wall shear stress, it 
was found that the shear stress decreases monotonically beginning from the 
entrance. Near the entrance, the drop in shear stress is also very rapid. 
It levels off as the flow moves downstream. Therefore, the rise in the 
friction coefficient is due to the rapid deceleration of fluid velocity 
downstream. 

The Nusselt number for plane flow under zero gravity also shows a 
rapid drop close to the entrance due to the development of the thermal 
boundary layer. After a short distance, however, the Nusselt number 
becomes nearly constant for the cases of heating and evaporation. As also 
seen by previous investigators, the Nusselt number for heating is larger 
than the Nusselt number for evaporation. A plot of the velocity vectors 
for plane flow under zero gravity is shown in Fig. 5.5. It is observed 
that the velocity profile is parabolic in nature and the magnitude of the 
velocity decreases downstream. 
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The results of radial flow under zero gravity is shown in Figs. 5.5 

and 5.6 for Re. =34. It can be seen that for the flow conditions 
in 

considered here the film height increases monotonically from the entrance. 
However, as demonstrated by Rahman et al. (1989a), the film height may 
increase or decrease in radial flow in contrast to plane flow where the 
film height always increases. The film height is basically a function of 
the inlet Reynolds number. The friction coefficient, in the case of radial 
flow, is found to decrease first, attain a minimum, and increase thereafter 
all the way to the exit. As discussed for the case of plane flow, the 
decrease of the friction coefficient near the entrance is associated with 
the development of the velocity profile from a uniform entrance condition. 
The rise of the friction coefficient after the minimum in this case is 
found to be much more rapid than that in plane flow. This is because the 
flow velocity decreases more rapidly due to the spreading of the fluid as 
it moves downstream. The Nusselt number for radial flow has basically the 
same kind of behavior as in plane flow. The values of the Nusselt number 
are different however, particularly in the case of evaporation where there 
is a loss of heat from the free surface. The vector plot in Fig. 5.5 shows 
that in radial flow the velocity decreases more rapidly downstream, while 
maintaining a parabolic structure across the film. It can be also noticed 
from Fig. 5.5 that the flow in a zero gravity environment is 

well- structured and no instability or separation is found to be present in 
the flow field. 

5.5.2 Film Flows Involving a Hydraulic Jump 

A hydraulic jump is found to be present under some flow conditions 
where a film flows along a horizontal or inclined surface under the 
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presence of a gravitational body force. Two different flow regimes are 
found to be present in this kind of flow: a supercritical flow before the 
jump, and a subcritical flow after the jump. Depending on the Reynolds 
number, the Froude number, and the length of the flow domain, the flow may 
be only supercritical, only subcritical or mixed in the computational 
domain. A mixed flow is the focus of our attention here, since a hydraulic 
jump can occur only in that kind of flow. Both plane and radial hydraulic 
jumps are considered in this section. The results for plane flow are 
demonstrated in Figs. 5.7 and 5.8. 

Figure 5.7 shows a plane hydraulic jump at Re = 355 and - ®-0. 
The Froude number at the outlet is assumed to be 1.0 to simulate a critical 
outflow condition. A critical outflow is present when a film moves off of 
the plate and experiences a free fall under the action of gravity. If the 
plate is extended beyond the computational domain, the outflow is simply 
subcritical and an appropriate Froude number needs to be specified. 

As seen in the figure, the film height increases slowly before the 
jump, experiences a rapid rise at the jump front and thereafter decreases 
gradually all the way to the exit. A wave- like structure is found to be 
present at the commencement of the subcritical regime. This shape was 
previously photographed by Francis and Minton (1984) for a hydraulic jump 
in a water channel, but could not be captured in previous numerical 
investigations. Looking at the velocity vector plot in Fig. 5.8, it can be 
noticed that the fluid separates from the plate due to the sudden rise of 
the free surface at the jump location. Most of the mainstream flow tries to 
orient itself along the free surface since a zero- shear condition exists on 
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Horizontal distance, m 

Fig. 5.7 Computed results for plane hydraulic jump (Re=355,Fr in =9) 



.rttfm&tw. 



z — Z2 

-> : S.0000E-01 n«/s. z 2 = 0 .0 15m (jump region) 



Fig. 5.8 Velocity vectors in plane hydraulic jump (Re — 355,Fr; n — 9) 
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that surface. This results in a recirculating flow near the solid wall. 
After a distance following the jump, the flow reattaches again to the 
surface and eventually develops into a parabolic profile. The flow is also 
seen to be parabolic in the supercritical regime before the jump. The 
separation of the flow from the solid wall also results in a negative shear 
stress near the jump location as seen in the plot of the friction 
coefficient in Fig. 5.6. The friction coefficient decreases rapidly near 
the entrance due to the development of the flow from a uniform entrance 
condition to a parabolic structure. The separation and reattachment points 
can be identified from this plot as the locations where the friction 
coefficient is equal to zero. After the flow reattaches, the friction 
coefficient increases slowly due to the reduction of the flow velocity due 
to friction as the fluid moves downstream. 

The plots of the Nusselt number are shown in Fig. 5.7 for the cases of 
simple heating and evaporation. In both cases, the solid wall was assumed 
to be at a uniform temperature higher than the ambient. In both cases, 'the 
Nusselt number decreases rapidly, comes to a minimum and then increases 
further downstream. 

A hydraulic jump during the radial spread of fluid along a horizontal 

plate is demonstrated in Figs. 5.9 and 5.10, where Re- =50 and Fr. = 

in in 

11.0. As in the case of a hydraulic jump in plane flow, a critical outflow 
condition is assumed. As seen in Fig. 5.9, the supercritical flow is 
sustained for a short distance, where the film height increases faster than 
the plane flow. At the jump, there is a rapid rise in fluid level, which 
is followed by a gradual decrease in the subcritical region. The velocity 
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Fig. 5.9 Computed results for radial hydraulic jump (Re; n =355,Fri n — 11) 





Fig. 5.10 Velocity vectors in radial hydraulic jump (Rejn — 355,Fri„ — 11) 
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of the fluid decreases downstream, due to the combined actions of radial 
spreading and friction from the solid wall. There is a short separation 
region as seen in both vector plots in Fig. 5.10 and the plot of friction 
coefficient in Fig. 5.9. The profile is approximately parabolic in both 
the supercritical flow before the jump and the subcritical flow following 
the reattachment. The rapid drop of the Nusselt number and the friction 
coefficient are associated with the rapid flow development near the 
entrance. As seen in plane flow, the Nusselt number reaches a minimum in 
the separated flow region. The Nusselt number rises slightly after the 
flow reattaches to the surface. 

5.6 CONCLUSIONS 

The results of a new methodology for the computation of the free 
surface flow of a thin liquid film is presented. The problems considered 
are plane and radial flows involving a hydraulic jump and plane and radial 

flows in a zero-gravity environment when no jump can be present. The new 

* 

method computes the whole flow field as a single domain problem 
irrespective of any regime change in the computational domain. This allows 
the determination of the details of the flow field and the friction and 
heat transfer behavior around a jump location. It was also found that the 
new method yields results that are more accurate than any other existing 
method for all of the problems considered here. 

In both plane and radial flows where a hydraulic jump is present, it 
was found that a distinct supercritical regime is present before the jump 
and a subcritical regime exists downstream of it. The very rapid rise in 
the fluid level at the jump location results in the separation of the flow 
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from the solid wall and the development of a backward circulation. The 
flow, however, reattachs to the surface after some downstream distance in 
the subcritical region. The back flow was found to be stronger in plane 
flow and at higher Reynolds numbers. The location of the separation and 
reattachment can be determined from the plot of the friction coefficient, 
which is negative in the recirculation region. The Nusselt number was 
found to have a minimum in the jump region. 

In the case of plane and radial flows in a zero gravity environment, 
it was found that the variation of the film thickness is monotonic and no 
separation of the flow is present. The friction coefficient decreases 
rapidly close to the entrance, attains a minimum and thereafter increases 
gradually downstream. The Nusselt number remains almost constant except 
for the entrance region, and is higher for heating than for evaporation. 
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Section VI 


ANALYSIS OF THE FLUID FLOY AND HEAT TRANSFER IN A TUN 
LiqUID Fill IN THE PRESENCE AND ABSENCE OF GRAVITY 

6.1 SUHARY 

The hydrodynamic and thermal behavior of a thin liquid film flowing 
over a solid horizontal surface is analyzed for both plane and radially 
spreading flows. The situations where the gravitational force is 
completely absent and where it is significant are analyzed separately 
and their practical relevance to a micro- gravity environment is 
discussed. In the presence of gravity, in addition to Reynolds number, 
the Froude number of the film is found to be an important parameter that 
determines the supercritical and subcritical flow regimes and any 
associated hydraulic jump. A closed- form solution is possible under 
some flow situations, whereas others require numerical integration of 
ordinary differential equations. The approximate analytical results are 
found to compare well with the available two-dimensional numerical 
solutions. 
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6.2. imODPCTIOW 


A thin film may be generated either by impingement of a jet on a 
solid wall or by discharge of fluid through a slot from a pressurized 
container. It may also occur during a melting or condensation process. 
Such situations may occur in the Space Shuttle for inflow to a 
propellant tank and in the absorber unit of a heat pump absorption 
system. The understanding of such flows in a micro- gravity environment 
is essential for the optimal design of fluid flow and heat exchange 
processes in a space vehicle. 

The inviscid flow of a liquid jet impinging on a solid wall is a 
classical hydrodynamics problem which is available in textbooks 
(Batchelor (1967)). It is concerned with irrotational, incompressible 
and inviscid flow, in which the effects of gravity and surface tension 
are neglected. One of the major attractions of this type of approach is 
that it can be handled using complex potential theory and therefore can 
be treated analytically. 

In nature, however, viscous effects become important, particularly 
when the thickness of the liquid layer becomes small. Also, gravity 
cannot be neglected entirely in most situations. Vatson (1964) analyzed 
the fluid mechanics of thin films produced by the impingement of a 
liquid jet on a flat horizontal surface under the action of gravity. By 
using the boundary layer approximations of the governing transport 
equations, analytical solutions using a similarity transformation along 
with the Pohlhausen integral method were derived. The analysis covered 
the regions where the boundary layer thickness is less than the film 
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height and where the film is totally engulfed by the boundary layer. 
The possibility of a hydraulic jump in such a flow was also anticipated. 
However, the analysis was applicable only to the supercritical flow 
before the jump. The height of the jump was predicted for any given 
location of the jump. The heat transfer counterpart of the impinging 
jet problem was considered by Chaudhury (1964). The energy equation was 
solved in closed fora including the effects of viscous dissipation by 
approximating the temperature profile with a fourth-order polynomial. 
Nusselt numbers for different values of the Prandtl number were 

presented. 

The impingement of a liquid jet in a gravity- free environment was 
presented by Labus and DeVitt (1978). They included the effects of 
surface tension, but entirely neglected the viscous forces. From a 
scaling analysis, it was shown that for a large jet Reynolds number, the 
flow can be approximated to be inviscid. Numerical solutions of the 
governing transport equations were obtained and compared with 
experimental measurements. Three distinct flow patterns of the jet were 
obtained which were classified in accordance with relative importance of 
inertia and surface tension. 

In the studies mentioned above, the investigators considered either 
inviscid flow or viscous supercritical flow up to the location of the 
jump. The fluid flow in the vicinity of the jump or in the subcritical 
region following the jump is also important from an engineering point of 
view. Moreover, all of the studies mentioned above are concerned with 
a thin film formed by impingement of a liquid jet. The driving 
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mechanism of these flows are flow rate and jet diameter which are 
characterized by jet Reynolds number. So, the findings in these studies 
cannot be applied readily to a situation where the film is formed by 
discharge from a pressurized container. This situation was considered 
in recent studies by Rahman et al. (1989a, 1989c, 1989d). A systematic 
numerical study of two-dimensional fluid flow and heat transfer in a 
thin liquid film in both plane and radially spreading flows was 
performed. The studies covered both zero and normal gravity 
environments. In the absence of gravity, no jump was found and the flow 
remained supercritical in the entire domain. In the presence of 
gravity, a jump was found under some flow conditions. In all of these 
studies, the flow field was computed numerically using a boundary- fitted 
coordinate system where the irregular free surface of the film was taken 
as one of the boundaries of the computation domain. The inertial, 
viscous and pressure forces were identified to be dominant and the 
surface tension was found to be negligible in most regions. The 
computational methodology presented by Rahman et al. (1989a and 1989d) 
(described in sections III and IV of this report) is termed the 
’pressure optimization method’ where the shape of the free surface was 
represented by an algebraic equation with two or more arbitrary 
constants. The constants were optimized using an exhaustive search 
which minimized the difference between the computed free surface and 
ambient pressures. The method was found to be robust and was extremely 
satisfactory for zero-gravity flows. However, when a jump was present 
in the computation domain, the pressure optimization method was unable 
to compute the jump as a single- domain problem. The subcritical and 
supercritical flows were computed separately and were matched at the 
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jump interface preserving the conservation of mass, momentum and energy. 


The two-dimensional numerical method presented by Rahman et al. 
(1989c) (described in section V of this report) is termed ’porous wall 
method’. In this method, the free surface was assumed to be a permeable 
wall where fluid particles could cross this boundary depending on the 
difference in fluid and ambient pressures. The shape of the surface was 
improved by successive iterations until the free surface conformed to a 
streamline where no penetration occurred. This method proved to be 
successful in handling the regions before and after the jump as a 
single- domain problem. The details of the flow structure in the 
vicinity of the jump were presented along with values of the skin 
friction and heat transfer coefficients. 

Even though the complete numerical solution for the flow field and 
heat transfer coefficient are available to us, the need for a simple 
systematic, analytic method to describe the flow is inevitable. The 
numerical results are limited to the flow parameters used for the 
computation, and cannot be extended in general. Moreover, an analytical 
expression is easier to use and implement in addition to the 
understanding of the limiting behavior of the flow. The present work is 
intended to analyze both zero-gravity flows and the flows where gravity 
is significant. The application of these results in a micro- gravity 
situation will be examined. Unlike previous analytical studies, it will 
cover both supercritical and subcritical regions and the jump connecting 
the two domains. In addition to analyzing the flow field, the heat 
transfer will be studied. 


209 



6.3. EQUATIONS OF IQTIOW 

The schematic of the problem under consideration is shown in Fig. 
6.1. A thin liquid film is flowing adjacent to solid heated wall. Two 
classes of flows are considered here. 

(a) Plane flow: where the film moves in a two-dimensional fashion 

along a plane horizontal wall. 

(b) Axisymmetric radial flow: where the liquid is introduced at the 

center of a circular horizontal plate and spreads out radially. 


The z(r) axis is directed along the longitudinal (radial) 
direction, and y-axis is directed normal to the plate. The velocity 
components in these two directions are w and v, respectively. The 
height of the free surface from the solid wall is denoted by 6 which 
varies with the longitudinal (radial) location of the plate. 

The equations governing the conservation of mass, momentum and 
energy for an incompressible constant property flow are given by 


DV _ 
Bt = 


i Vp + i/V 2 v + g 


( 6 . 2 ) 


- a V 2 T 


(6.3) 


Here we have also assumed that heat generation due to viscous 
dissipation and pressure work are negligible, and there is no generation 
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of heat or mass. The velocity vector V can be resolved into components 
as 


V = vk + vj 


The vectorial representation used in equations (6. 1-6.3) is easier to 
handle since it can be readily transformed to Cartesian or cylindrical 
components as needed. The boundary conditions in component form are 
given by 


II 

< 

for isothermal wall 

(6.4) 

•sTs 

II 

< 

for constant flux wall 


d£ v 

= w’ p 

* Pe> T e = 0 


| T = T sat’ 

for evaporation 

(6.5) 

q e = 

for adiabatic condition 



at z = 0 or r = r. : 
in 



( 6 . 6 ) 


at z = L or r = r^: 



fa _ 0 
TR - °» 

Fr = Fr 



for zero gravity 
, , for non- zero gravity 


(6.7) 


On the free surface, both streamline and stress- free conditions 
have to be satisfied. The balance of normal stresses, in general, 



relate fluid pressure to the ambient pressure via surface tension and 
other stresses. From a scaling analysis, it can be shown that (see 
Rahman et al. (1989c)), for a reasonably large Veber number and flow 
rate that is typical for these flows, these stress terms are found to be 
an order of magnitude lower than the pressure. So, the surface tension 
can be assumed to be negligible in most regions of the flow leading to 
the p = p e condition on the free surface. 

Both plane and radial flows described above will be solved for zero 
and normal gravity situations. In the following two sections, we will 
analyze the flow using uniform and parabolic velocity distributions, 
respectively. The former is more common in fluid flow literature 
concerning a hydraulic jump or shock wave and will be carried out in 
details. The latter section will improve on the analysis by using a 
parabolic velocity distribution across the film, which is more 
appropriate for laminar thin film flow. 

6.4 ANALYSIS OF FLOY USING ONE- DHENSIONAL UNIFORM VELOCITY 

Ve first consider the situation where the velocity variation across 
the thickness of the film is neglected, and the film is assumed to 
propagate downstream with its uniform average velocity. Let V be the 
average velocity of the film in the longitudinal (radial) direction of 
the plate and Q be the volumetric flow rate. In the case of plane flow, 
Q is the volume flow rate per unit width. The continuity equation (6.1) 
can be written as 

q = (2 t r) A SV (6.8) 
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, , f 0, for plane system 

where ^ = | 1} for rad ial system 


Integrating eqn. (6.2) across the thickness of the film, expressing 
the resistance from the solid wall in terms of friction coefficient and 
substituting eqn. (6.8) results in 


3F Cr * sO 


Cf V 2 

■ r r 


(6.9) 


This momentum equation must be solved along with eqn. (6.8) to 
determine the flow field, and will be carried out in the following 
subsections . 


6.4.1 Flow Under Zero Gravity 

For a steady flow under zero gravity, the governing equations (6.8 
and 6.9) reduce to 

Q = (2rr)^ <5V = constant (6.10) 


dV C f V 2 

3? = ' T T 


Eliminating 6 from eqns. (6.10 and 6.11) gives 


( 6 . 11 ) 


V 2 dV - - (2 xt) dr 

Integrating this equation assuming a constant c^ and substituting 
the conditions at one location of the flow (i.e., at r = r in , V = V ifl 
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f 


and 6=5-) yields the solution in the form 
in' 


V 

VT 

m 


t (A +~ry i 


1 - 


' in 




+ T) <5^ 


5 in / r in A 
= h ' r ' 


m 


Solving for 6 results in the following relation 


« = i ^ * (i . i> ( !ii) A 

in in 


Cr r . 

f in 


where A - T ^ 


( 6 . 12 ) 


(6.13) 


(6.14) 


{ From equation (6.13) we may express the distribution of film height for 

r 

plane and radial flows. 

i 

! 

For plane flow (A = 0), the equation simplifies to 


( 


i 



1 + 


c f 

T~ 



(6.15) 


This indicates that for a constant friction coefficient the film height 
increases linearly with distance. This increase in film height is 
because of the decrease in flow velocity due to the resistance from the 
wall. Ve may also notice that for inviscid flow when no resistance is 
exerted by the solid wall, the film height remains the same at all 
downstream locations beginning from the entrance. 
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In real flow situations, however, the friction .coefficient varies 

i 

along the plate. Usually, near the entrance the friction coefficient is 
greater due to the rapid changes in the velocity profile as the boundary 
layer develops beginning from the entrance plane. A good estimate of 
the friction coefficient can be made using the Blasius solution (see 
Kays and Crawford (1980)) for forced convective boundary layer adjacent 
to a plate. This results in 


c 


f ' 


0.664 

: 7z~ 

v 

* 


for plane flow 


j _ 0.664 

and C| = ■ — , 


V (r-r. T 
v in' 


for radial flow 


(6.16a) 


(6.16b) 


However, unlike any other external flow, in a thin film the boundary 
layer thickness is frequently of the order of the film height. So, 
after some distance downstream, the viscous effects propagate all the 
way to the free surface. Then can be better estimated from the 
equation presented in a later section of this paper (equation 6.46). 
The location where one should switch over from one formulation to the 
other may be estimated from the solution of the boundary layer thickness 
for forced convection. The thickness of the momentum boundary layer is 
given by 
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5.0 



vz 

Jr 


, for plane flow 


and 



for radial flow 


As shown by Thomas et al. (1989), the two friction factor estimates 
yield the same value at S = 1.81 So, for continuity, one may use 
equation (6.16) for 6/6^ > 1.81 and equation (6.46) for 6/6^ < 1.81. 


The one- dimensional analytical solution with the friction 
coefficient described above has 'been compared with the numerical 
solution of Rahman et al. (1989a) in Fig. 6.2 for Re = 12.5, L = 0.03 m 
and <^ n = 0.000595 m for plane flow. Similar to the analytical 
solution, a uniform inlet velocity has been assumed for the numerical 
solution. The predicted film height distribution is found to be 
somewhat higher than the two-dimensional numerical solution. A test run 
with the friction coefficient reported in Rahman et al. (1989a) yielded 
a solution almost coincident with the numerical film height. So, the 
discrepancy is primarily due to approximations inherent with the 
estimation of the friction coefficient. The Blasius solution is 
perfectly valid for boundary layer flow where the free stream extends 
through a large distance and no free surface is present. However, the 
comparison is reasonable considering the approximations inherent in the 
simple analytical formulation. 


A plot of equation (6.13) for radial flow under zero gravity is 
presented in Fig. 6.3 for different values of A. For A = 0.0, which 
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Fig. 6.2 Film height for plane and radial flows under zero gravity 




corresponds to inviscid flow, the film height monotonically decreases 
with radial distance. This is quite expected since, in the absence of 
any friction, the velocity of fluid particles will remain constant. 
Therefore, the film height has to decrease with radius as more and more 
area becomes available to the flow. In this plot we also notice that 
when friction is present (A > 0), the film height may increase or 
decrease depending on the value of A and the radial location. Also, for 
a given radial location the value of the film height is more for a 
larger value of A, since A quantifies the frictional resistance compared 
to the inertial forces. For A = 0.2 and 0.4, the film height decreases 
first, attains a minimum and then increases further downstream. The 
location for minimum 8 can be determined by differentiating equation 
(6.13) which gives 

(6.17) 

(6.18) 

From eqn. (6.17) we find that a minimum exists if A < 1/2. Otherwise, 
the film height increases continuously from the entrance as seen in 
plots corresponding to A > 0.5. At A = 1.0, the film height increases 
linearly with radius. 

The analytical solution for a particular case of Re in = 404, $ in = 
0.005 m, r. = 0.05 m and r . = 0.2 m is shown in Fig. 6.2 where it is 
compared to the two-dimensional numerical solution given by Rahman et 
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al. (1989a). Here, we have also used the friction coefficient estimated 
by equations (6.16) and (6.46). The results are very close except for 
the fact that the analytical solution shows a minimum in the flow 
domain. The minimum is indeed present as verified from the value of A. 
The numerical solution could not accommodate this behavior since a 
continuous hyperbolic curve was assumed to represent the free surface. 

6.4.2 Flow in the presence of gravity 

Vhen the gravitational body force is significant compared to other 
existing forces, its effect should be included in the analysis. Unlike 
a regular forced or natural convection flow in outer space, the effect 
of gravity, even small, may be significant in a thin film flow since two 
potential flow regimes, namely supercritical and subcritical may be 
encountered. These regimes are characterized by the Froude number. The 
transition of the flow from supercritical to subcritical is analogous to 
the transition from supersonic to subsonic flow in gas dynamics, where 
the Mach number determines which regime is present. The analysis 
presented below uses the methodology commonly followed in the analysis 
of high speed flow with friction (i.e., Fanno flow). The application of 
these ideas to a thin film flow is entirely new and it enabled us to 
combine both flow regimes (supercritical and subcritical) in a single 
analysis which was not possible by previous investigators. Most 
previous analytical studies assumed the film height to remain constant 
in the subcritical region. Expressing eqns. (6.8) and (6.9) in terms of 
the Froude number results in 

q = {2tt) X /g Fr tf 3 / 2 ( 6 - 19 ) 
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( 6 . 20 ) 


Sr [* ( 1 * T » * ' ^ Ff2 

These equations will be solved for plane and radial flows in the 
following subsections. 

(a) Plane Flow 

For plane flow ; = 0. Then eqns. (6.19) and (6.20) can be combined 
by eliminating the film height and non- dimens ionalized assuming c f to be 
locally uniform to give the equation 

(1 ■ Fr 2 ) jj|f = | Fr 11 / 3 (6-21) 


where 



7^773 


( 6 . 22 ) 


Integration of equation (6.21) gives 


^ Fr' 8 / 3 - j Fr' 2 / 3 = - R1 + c 

To evaluate the integration constant, the Froude number must be 
specified at one location in the flow. Upon examining equation (6.21) 
we notice that the equation is singular at Fr = 1, so a critical flow 
situation is present at that location. This is analogous to a Mach 
number of unity in a compressible flow. Let this critical location be 
denoted by Rl*. Then the solution can be written as 
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(6.23) 


u-u*. - I Fr-*/* - J, ft-*/* 

Equation (6.23) is a double- valued function as shown (in circles) 
in Fig. 6.4. The two branches of the function represent subcritical or 
supercritical flows where the Froude number is less than or greater than 
unity, respectively. 

Since two solutions exist at any location, the possibility of a 
sudden jump from supercritical to subcritical flow exists. The opposite 
is not true since that would violate the second law of thermodynamics. 
The height of the film before and after the jump can be related by the 
conservation of mass and momentum across the jump. This is given by 

^ = j [ J 1 * 8 Fr^ - 1 ] (6-24) 

where subscript 1 indicates conditions before the jump and subscript 2 
indicates conditions after the jump. 

In Fig. 6.4, it can be observed that both subcritical and 
supercritical flows move towards the critical condition. In the 
supercritical regime, the Froude number decreases downstream, whereas in 
the subcritical regime it increases with distance. Since the Froude 
number is inversely related to the film height, the film height is 
expected to increase downstream in the supercritical region and decrease 
downstream in the subcritical region. A flow starting with Frcl will 
follow the subcritical curve all the way to the exit. However, a flow 
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Fig. 6.4 Analytical solution for plane and radial flows over 
a plate in the presence of gravity (uniform velocity) 


starting with Fr>l may follow the supercritical curve all the way to the 
exit or may encounter a jump and transform to subcritical flow inside 
the region under consideration. The location of the jump is determined 
by the downstream condition of the flow. 


(b) Axisvmmetric Radial Flow 

For radial flow, A = 1, eqns. (6.16 and 6.17) may be transformed 
into the following form 


(1 


t>_2\ dFr 

Fr > 3ET 


9 „_ll/3 Fr(2 + Fr 2 ) 

2 Fr 2 Ri 


(6.25) 


where 


Rl 


r 


(i/ 2 / g ) 1 / 3 ae 2 ^ 3 

C f 


(6.26) 


Analogous to plane flow, here the radial distance has been 
non- dimensionalized using the local Reynolds number and friction 
coefficient. The variation of friction coefficient with radial distance 
has been assumed to be negligible. A closed- form solution of equation 
(6.25) is not possible, so a numerical integration was performed using 
the Euler method. It can be noticed that the equation is singular at Fr 
= 1. Therefore, the critical condition cannot be directly applied as a 
boundary condition in the numerical solution. To avoid this 
singularity, the equation can be expanded around the singular point and 
the solution can be found at a short distance from the singular point 
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from the lowest- order expansion. The numerical integration then can be 
carried out beginning from a short distance away from the singular 
point, where the solution is already known. The solution is shown in 
Fig. 6.4. 

♦ 

It should be noted that the critical radius, R1 , appears as a 

parameter. For a given flow rate and inlet Froude number, the value of 

* * 

R1 depends on the inlet radius. A large R1 indicates that the flow 

starts up at a larger radial distance, where the effect of curvature is 

* 

smaller. This is quite evident from the plot since the result for R1 = 
100 coincides with the solution for the plane flow. The double- valued 
nature of the solution is also present in radial flow which also 
indicates the possibility of a hydraulic jump. The film height before 
and after the radial jump can be related by the same equation as for 
plane flow (equation 6.24). 


6.4.3 Characteristic Behavior of the Flow 

Since the equations of transport for the free surface flow of a 
thin liquid film are somewhat similar to those for one- dimensional 
compressible flow, it may be useful to analyze the characteristic 
behavior of the flow. 


The conservation equations in its time- dependent form can be 
written as 


1 3 
m + 3 ? 



E) 


= H 
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where 



c 



u = 

0 

6V 

, E = 

S(V 2 ♦ §*) 


and H = 


V' 


where r y is the surface shear stress. These are two first-order partial 
differential equations in t and r with two dependent variables, 5 and V. 


These two equations may be linearized and written in the following 
characteristic form: 


R t . C R x = S 


where 


V * J- 


= Riemann invariants 


C = 


* /glj 

- r&i 


- Vave speed 


It can be seen that the first invariant always propagates downstream 
(i.e., C > 0) . The second invariant, however, propagates downstream for 
supercritical flow (Fr = V^//~g7^ > 1) and propagates upstream for 
subcritical flow (Fr < 1) . This implies that both V and 6 must be 
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prescribed upstream for solving supercritical flows while only one is 
prescribed upstream and one downstream for solving subcritical flows. 


6 - 4 - 4 Application of Results to Mir.ro- Gravity 

The zero-gravity thin film flow considered here can be attained 
only in a complete absence of the gravitational body force. However, in 
orbit the gravity is very small, but not precisely zero. Therefore, the 
range of applicability of the results needs to be investigated. Looking 
at plane flow results under zero gravity we notice that 

S = ( tS- ) z + 8- 
K 2 J in 

This suggests a monotonic increase in the film height so long as c f 
remains positive. The flow decelerates due to friction and consequently 
the film height increases. Since the major driving mechanism for this 
flow is inertia, a film introduced with a finite velocity becomes very 
slow after traveling a certain distance. 

Vhen investigating the expression for the Froude number, it can be 
noticed that even for a very small gravity force, where the inlet Froude 
number tends to infinity, the Froude number may become on the order of 
unity after some distance, since both the deceleration of flow and the 
increase of film height contribute to reducing the Froude number. The 
situation here is analogous to hypersonic flow, which in the presence of 
friction rapidly ceases to be hypersonic after traveling some distance. 
Therefore, in a microgravity situation, the order of magnitude of the 
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local Froude number is very important. Vhen the Froude number is 
extremely large, like near the entrance, the zero gravity results may be 
applicable. However, away from the entrance, one should resort to the 
normal gravity results presented here and eventually a hydraulic jump 
must be encountered. 

6 # 5 ANALYSIS OF FLQV AND HEAT THAN SF El USING PARABOLIC — VELOCITY 

DISTRIBUTION 

The one- dimensional analysis presented in the previous section 
revealed many interesting features of the flow. However, the velocity 
profile in reality is two-dimensional in nature because of the no- slip 
condition at the solid wall. As demonstrated by Rahman et al. (1989a), 
the velocity profile is approximately parabolic in nature in most 
regions of the flow. The temperature profile may also become parabolic 
after the thermal boundary layer develops. 

Ve assume the velocity component, w, and temperature, T, to have 
the following general form where the boundary condition at the solid 
wall and free surface are satisfied. 

}= 3 ^ (1 - ( 6 - 27 ) 

( 6 - 28 ) 

V 

where J] is the dimensionless coordinate across the thickness of the 
film. The constants B and C will be evaluated later in this section for 
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different thermal conditions considered here. 


Once the velocity and temperature profiles are fixed, the other 
flow quantities can be readily calculated and are given as follows 

r v = ( 6 - 29 ) 

KBT 

9 w - ~~S (6.30) 

KBT 

9 e = ' ~y W (® + 2 ^) (6.31) 


T e = T w (i + B + C) (6.32) 

T b = T w + 5 /8 B + 9/20 C) (6.33) 


In a thin film flow, the velocity across the thickness of the film 
is much smaller than the velocity along the plate. Scaling the governing 
equations (6. 1-6.3) using the condition v << w gives 


1 d 
^ 2F 


< * ^ 
(r w) 



(6.34) 


w 


dv 

3F 


+ 


dv _ 1 dp 

V ~yr-~ — - — -rr*- + 

(7Y p or 


02 
a w 

v — 5 


( 6 . 35 ) 


i & + 

p o y 


g - 0 


( 6 . 36 ) 
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(6.37) 


an an _ a 2 i 

Here A = 0 indicates plane flow and A = 1 denotes radial flow. 

Integrating equation (6.36) we get, 

P = P e + Pg (<5-y) (6.38) 

where p g is the pressure at the edge of the free surface. This 
indicates that pressure across the thickness of the film is hydrostatic 
in nature. Also, at zero gravity, p = p g . In the external flow over a 
flat surface, the ambient pressure usually remains constant. So, in a 
gravity- free environment, the static pressure is expected to remain 
constant everywhere in the flow field. 


6.5.1 Plane Flow Under Zero Gravity 

For plane flow under zero gravity equations (6.34) and (6.35) can 
be integrated from 0 to 6 using the velocity profile given in equation 
(6.27) and the no- slip condition at the solid wall and zero- shear 
condition at the free surface. After using the definition of Reynolds 
number, it results in the relationship 


d <5 _ 2.5 

= w 


(6.39) 


Integrating this equation, one obtains 


6 

T~ 


in 


1 + 


ir 



(6.40) 
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For plane flow, the Reynolds number remains constant. Therefore, 
the film height increases linearly beginning from the entrance location. 
Ve also notice that in the limiting case of Re -* «, 5 = ^ in everywhere 
in the flow field. This is quite expected since in an inviscid plane 
flow, the film height remains constant. 

Figure 6.2 shows a comparison of this solution with previous 
numerical results of Rahman et al. (1989a) . The analytical solution 
tends to predict a somewhat higher film height than the numerical 
solution, but lower than the solution predicted by assuming a uniform 
one- dimensional velocity distribution. The assumption of a parabolic 
velocity profile, even though not exact everywhere in the flow field, 
produces a reasonable estimate of the film height distribution. 

6.5.2 Radial Flow Under Zero Gravity 

For radial flow under zero gravity, the integrated equation for the 

conservation of momentum can be written as 

d 6 2.5 6 (6.41) 

3r = Ee~ ' r 

Analogous to the plane flow case, we have used the continuity 


equat 

ion, the boundary 

cond 

itions at the 
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and 
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resistance at the solid wall reduces the flow velocity and consequently 
tends to increase the height. However, at the same time, the flow is 
spreading radially with more area available for the flow as the radius 
increases. This same behavior was seen in the discussion of 
one- dimensional analysis (Fig. 6.3), where the parameter ’A’ determined 
the variation of film height downstream from the entrance. A is a 
function of inlet Reynolds number along with inlet height and radius. 
Note also that the Reynolds number does not remain constant in radial 
flows and changes with location. 

Integrating eqn. (6.41) gives 


5 *i. 2.5 (r 3 - 

3 OH “in r in 


( 6 . 42 ) 


The location of the minimum film height may be calculated by 
differentiating equation (6.42) with respect to r. This results in 


r 3Re- h- 


( r )for min 6 r in 


in“in / 1 5 ^ 

5— ( f7: ' BUe-TT-J J 


m 


1/3 


( 6 . 43 ) 


Ve can also notice that the minimum will exist when 


Re • h • r 
in in , 5 

< j 


m 


The results of the analytical solut ion for radial f 1 >>v unde: zero 
gravity is also compared with previous numerical results of Rahman et 
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al. (1989a) in Fig. 6.2. In this case, we can see a minimum inside the 
flow domain confirming the observation from the one- dimensional 
analytical solution using the uniform the velocity profile. Both 
analytical solutions are close to each other whereas the numerical 
solution is slightly higher in most parts of the flow. 


6.5.3 Plane Flow in the Pr esence of Gravity 

For plane flow where the gravitational body force term is retained, 
the governing equations of motion (6.34-6.35), can be integrated across 
the thickness of the film using equation (6.38), the boundary 
conditions, and the definition of Reynolds number, and can be 
transformed to the equation 



where 



r 2,1/3 


g J 


Re 


5/3 


(6.44) 


(6.45) 


Here the radial coordinate R is normalized in terms of Reynolds number. 
Note that the definition of R is very similar to the normalized radius 


R1 used in the one- dimensional analysis using the 
except that the friction coefficient c^ appears in 
not appear in R. From equation (6.29), the friction 
easily calculated to be 


uniform velocity 
R1 whereas it does 
coefficient can be 
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(6.46) 


This friction coefficient is correct when the velocity profile is 
parabolic across the thickness of the film, which obviously assumes that 
viscous effects are extended all the way to the free surface. Vhen this 
definition of C£ is assumed, III becomes the same as R. In the present 
investigation we have tried to keep the one- dimensional analysis using 
the uniform velocity somewhat more general where any known distribution 
of the friction coefficient may be used. 

From equation (6.44) it can be noticed that a critical condition in 

the flow occurs at Fr = 0.913. Note that this value of Froude number is 

different from the conventional critical Froude number of Fr = 1, which 

is strictly valid for a one- dimensional flat velocity distribution. The 

value of Froude number of 0.913 is therefore the "weighted" value for 

the parabolic profile. Integrating equation (6.44) with a boundary 

* 

condition of (R = R at Fr = 0.913) we can obtain the solution as 

- 1.195 + I Fr' 2 / 3 - ^ Fr' 8 / 3 = 3.75 (R - R*) (6.47) 

A plot of this equation is shown in Fig. 6.5 (as circles). Analogous to 
the one- dimensional analysis, the solution is also found here to be 
double valued in nature. The two stems of the curve denote the 
supercritical and subcritical flow regimes. The possibility of a 
hydraulic jump also exists here, which basically depends on the incoming 
froude number, length of the plate and flow condition at the outlet. 

The jump may be present only when the flow initially is supercritical. 
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J-l 


Fig. 6.5 Analytical solution for horizontal flow over a plate 
in the presence of gravity (parabolic velocity) 



The ratio of the film height before and after the jump can be determined 
from a mass and momentum balance at the jump, which, in this case, turns 
out to be 




1 

5 


1 + 9.6 Fr^ 


1 


(6.48) 


where the subscript 1 indicates conditions before the jump and the 
subscript 2 indicates conditions after the jump. In comparing this 
relationship with equation (6.24) for uniform flow only the coefficient 
within the radical changes for the parabolic profile. 


Plane flow in the presence of gravity is characterized by two 
independent dimensionless groups, namely, the Froude number and Reynolds 
number. Since the Reynolds number is lumped with other parameters in 
the definition of R, to see its effects more clearly, the plane flow 
solution is plotted again in Fig. 6.6 for different values of Reynolds 
number using (z - z )/(i^/g) ^ as the abscissa. Note that for a small 
value of Reynolds number, the supercritical solution cannot be sustained 
for a large length of the plate. On the other hand, for the same 
Reynolds number, the subcritical solution may be present for any length 
of the plate. Also, the length for which a supercritical solution may 
be present increases with Reynolds number. So, for a given length of 
the plate and inlet Froude number, the length of the supercritical flow 
regime is directly dependent on the Reynolds number. lor a ..iige 
Reynolds number, the flow may be supercritical in the onl i:o d"Oiin. 


Vit.h a decrease of the Reynolds number, a jump is expec 


,ed 


t o 


iear 
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F'ig. 6.6 Froude number distribution for plane flow over a plate 

in the presence of gravity 


with a subcritical regime. As the Reynolds, number is reduced further, 
the jump may move all the way to the inlet, resulting in completely 
subcritical flow. 

The figure also shows that if the flow enters the control volume 
with a higher Froude number, it may remain supercritical for a larger 
distance. The effect of gravity can be also analyzed. For a given flow 
rate and Froude number, (z - z ) is inversely proportional to (g)^/^. 
So, for a smaller value of g, the length of the supercritical flow 
regime is expected to be longer. In the limiting case of zero gravity, 
a subcritical flow regime is not possible and the flow becomes entirely 
supercritical . 

The predicted film height is compared with the previous numerical 
solution of Rahman et al. (1989a) for Re = 89, Fr. = 8.6, Fr A = 1.0 
and L = 0.14 m in Fig. 6.7. The height in the supercritical region is 
slightly over- predicted, whereas in the subcritical region, the height 
becomes flatter than the numerical solution. In the numerical 
computation of Rahman et al. (1989a), a critical outflow condition is 
assumed to be present at the exit, whereas in the analytical solution 
this condition was not imposed. 

6.5.4 Radial Flow in the Presence of Gravity 

For radial flow in the presence of gravity, the equations of motion 
(6.34-6.36) can be expressed as, 
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SB 


(6.49) 



dFr 5Fr + 3Fr 3 , 15 P „ll/3 

air = si — r ¥v 


where R = 



(6.50) 


Here, we can also notice that a critical condition is arrived at Fr 

= 0.913. This confirms that a critical condition in the flow depends on 

the velocity profile used for the analysis. The double- valued nature of 

the solution and the possibility of a hydraulic jump still exists. The 

* * 

solution is presented in Fig. 6.5 f<?r different values of R . Analogous 

to plane flow, the critical radius R corresponds to Fr = 0.913. As 

discussed before, a smaller value of R indicates that the flow starts 

up at a smaller radius where the effects of curvature are important. 

* 

This fact is quite evident in the figure. The curve for R = 100 
coincides with the plane flow solution. At this situation, the effects 
of curvature is negligible. Also, at R* = 0.5 and 1, the Froude number 
attains a minimum in the subcritical regime. This indicates that if a 
jump happens before that radial location, the film height may still 
increase in the subcritical region, attain a maximum and then diminish 
further downstream. This phenomenon is not present in plane flow and at 
flows with large values of R , where the film height decreases 
monotonically in the subcritical region. It can also be noticed that 
equation (6.49) reduces to equation (6.44) as R -• ®. So, the plane flow 
may be treated as a limiting case of radial flow where R -• a>. 


6.5.5 Analysis o f Heat Transfer 

The heat transfer behavior can be solved by integrating equation 
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of equation (6.28) for the temperature 


(6.37) with the substitution 
profile . 


Different thermal boundary conditions are considered in the present 
study. They are isothermal and uniform flux conditions at the solid 
wall and evaporation and pure heating without evaporation on the free 
surface. The heat transfer coefficient for these cases are defined as 


h = 


q w / (T w - T b ), for heating 


(6.51) 


[% / (T - ^sat^’ ^ or eva P orat ^ on 

A general expression for the Nusselt number can be evaluated from 
equations (6.30-6.33). 


* 


Nu = 


B 

q , for heating 

B 

B + C , for evaporation 


(6.52) 


A complete analytical solution is possible for the case of heating, 
when the free surface is assumed to be adiabatic in nature. For this 
case, Nu* = 2.5 for both isothermal and uniform flux wall conditions. 


* 

Note that Nu is a 
height, both of which 
oi this result with the 
( 1989d) is shown in 


function of heat transfer coefficient and film 
change with location on the plate. A comparison 
previous numerical solution of Rahman et al . 
Fig. 6.8. For radial flow at zero gravity (Re = 




'V-'i J 


z /L or (r r j n )/( r out r in) 

Fig. 6.8 Heat transfer coefficient for heating with no 
evaporation adjacent to a constant flux wall 



404, Pr = 7, r in = 0.05 m, r Qut = 0.2 tn, <S in = 0.005. m, T in = 20°C, q y = 

1000 V/m 2 ), the numerical Nusselt number starts at a higher value and 

then approaches an asymptotic limit as the flow moves downstream. This 

limiting value is somewhat lower that the analytical prediction. This 

indicates that the true temperature profile is somewhat flatter than the 

parabola assumed here. For plane flow in the presence of gravity (Re 

89, Fr- = 8.6, Fr . = 1.0, Pr = 7, L = 0.14 m, T in = 20°C, q y = 1000 

in oux 

V/m 2 ), the subcritical and supercritical heat transfer coefficients 
approach the analytically predicted value as the flow moves downstream 
with the development of the temperature profile. 


For an isothermally heated wall with evaporation on the free 
surface, an expression for the Nusselt number can be derived by 
integrating equation (6.37) using equation (6.28) for the temperature 
profile and boundary conditions. The result is 

»„* = !- — i. — e^*) (6-53) 

|sat .1 


where 


P( 2 ) 


80 a f Z dz 

nl,r 


(6.54) 


Here Cq is the value of C at one location in the flow, where integration 
should start. 


gravity 


expressions for Nu are true both for the cases where the 

is finite and infinitesimally small. The actual heat transfer 
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coefficient, however, will be different^ because of different 
distributions in the film height. 


6.6 conclusions 

A systematic analysis of thin film flows in a zero gravity and 
non- zero gravity environment is performed and their implications in a 
real flow problem under micro- gravity conditions is discussed. The 
analysis is carried out in two parts. First, the hydrodynamics of the 
flow is studied for uniform one- dimensional velocity with any given 
friction factor. Next, the flow field and heat transfer coefficient are 
determined using the Pohlhausen integral method. 

It is found that in the complete absence of gravity, the flow is 
supercritical, where even for a relatively small, but not negligible 
gravity, two potential flow regimes - supercritical and subcritical may 
be present. The two regimes are separated by a jump where large changes 
of film height take place. In the presence of friction, both 
supercritical and subcritical flows move towards a critical condition. 
The flow regime is characterized by the Froude number, whereas the 
Reynolds number indicates the relative importance of the inertial and 
viscous effects. 


For plane flow, it is 
film height increases linearly 
the film height increases 
encounters a jump and then 
subcritical region. The heat 


found that in the absence of gravity, the 
with distance. When gravity is present, 
monoton i rally in the supercritical region, 
decreases furtner downs '.ream in the 
transfer coefficient in this situation is 
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characterized by Nu* = 2.5, both for isothermal and uniformly heated 
walls with no evaporation from the free surface. 

For radial flow it is found that in addition to the Froude and 
Reynolds numbers, the inlet radius is an important parameter, which 
carries the effects of curvature in the flow. At large values of the 
inlet radius, a plane flow situation is approached. Here, in the 
presence of gravity, the film height may increase or decrease depending 
on the flow rate, radial location and Froude number. The comparison of 
the results with previous numerical solutions shows reasonably good 
promise that the analytical method presented here can act as a useful 
tool to easily obtain approximate results for any given flow situation 
in a normal gravity or microgravity environment. 
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Section VII 


COMPARISON OF THEORETICAL PREDICTION AND EXPERIMENTAL DATA 


7.1 SUMMARY 

The theoretical models developed in earlier sections are improved and 
a systematic procedure to compute the free surface flow of a thin liquid 
film is suggested. The solutions for axisymmetric radial flov on a 
stationary horizontal disk and that when the disk is rotating around its 
axis are presented. The theoretical predictions are compared vith the 
experimental data presented in Section II of this report. The details of 
the flov field are also presented for the stationary disk, which shows flow 

separation near the location of the jump. The effects of surface tension 

are found to be important near the outer edge of the disk where the fluid 


experiences a free fall. At other locations, the surface tension is 
negligible. For a rotating disk, the frictional resistance in the angular 
direction is found to be as important as that in the radial direction. 
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7.2 INTRODUCTION 


The appropriate design of a rotating absorber unit for a heat pump, a 
rotating condenser or evaporator, or the proper method of spin coating 
deposition on a metal surface all require an understanding of thin film 
flows over a horizontal rotating surface. Moreover, the impingement of a 
liquid jet on a solid wall that leads to the formation of a thin film is 
found in many engineering devices. These free surface flows are sometimes 
associated with a hydraulic jump where two different flow regimes are 
encountered at the upstream and downstream sides of the jump. The flow 
structure around a hydraulic jump is usually complicated in nature, the 
understanding of which is essential from both theoretical and practical 

points of view. 

The radial spreading of a liquid jet impinging on a flat horizontal 
surface was studied by Vatson (1964). Both analytical and experimental 
studies of laminar and turbulent flows were made. A hydraulic jump was 
found to be present under some flow conditions. The analysis covered the 
supercritical region before the jump where four different flow regimes were 
identified and analyzed using the Pohlhausen integral method and a 
similarity transformation. In the subcritical region, the film height was 
assumed to be constant. For a given location of the jump, the subcritical 
height could be also predicted. 

The heat transfer to a thin film formed by an impinging liquid jet was 
considered by Chaudhury (1964). His analysis basically followed the ideas 
introduced by Vatson (1964) except that the energy equation was solved. In 
the region away from the impingement location where a similarity solution 
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of the momentum equation was possible, the energy equation was solved in a 
closed form including the effects of viscous dissipation, in approximate 
solution for the simultaneously developing hydrodynamic and thermal 
boundary layers were also obtained using the integral method where the 
temperature profile was approximated by a fourth- order polynomial. Nusselt 
numbers for different values of the Prandtl number were presented for both 

the regions. 

The laminar flow of a thin film adjacent to a horizontal rotating 
surface was studied b, Rauscher et al. (1973). in asymptotic expansion 
technique was used to develop an approximate solution where the radial 
spread of fluid was perturbed to determine the effects of convection, 
Coriolis acceleration, radial diffusion, surface curvature and surface 
tension. The physical significance of these higher-order effects were also 

discussed . 

The transient behavior of the film when the rotation is started from 
an initial stationary condition was examined by Higgins (1986) . Both low 
and high Reynolds number flows were considered. An asymptotic expansion 
technique was used to develop analytical solutions valid for short and long 
time scales. Bis results were valid only when the thickness of the 
boundary layer is small compared to that of the film. The thin film 
of a non- Newtonian fluid adjacent to a rotating surface was considered in 
the study of Matsumoto et al. (1982). An integral method was used where 
the velocity profile across the thickness of the film was approximated by a 

polynomial . 
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The condensation of a saturated vapor onto a rotating horizontal 
surface was investigated by Sparrow and Gregg (1959). The governing 
momentum and energy equations were simplified to a set of ordinary 
differential equations by using a similarity transformation. The equations 
were then integrated numerically to predict the condensate layer thickness 
and heat transfer coefficient along with the temperature and velocity 
profiles. The evaporation of a thin film on a rotating surface was studied 
by Butuzov and Rifert (1973). Closed-form solutions were derived for the 
film thickness and heat transfer coefficient. 


In all the studies mentioned above, the investigators have developed 
analytical models for the phenomenon of the radial spread of fluid with or 
without rotation. Most of these models are for a thin film generated by 
the impingement of a liquid jet which is characterized by the jet velocity 
and its diameter. Therefore, they cannot be readily applied to a situation 
where the film is generated by a discharge from a pressurized container in 
the center of the disk. Moreover, none of these models handle very well 
the hydraulic jump associated with these flows. The studies by Vatson 
(1964) and Chaudhury (1964) analyze in detail the supercritical region 
before the jump, but make a simple assumption that the film height after 
the jump remains constant, which is not true in reality. A systematic 
investigation of thin film flows discharged from a pressurized container in 
the center of the disk has been carried out recently by Thomas et al. 
(1989) (described in section I of this report), Rahman et al. (1989a) 
(described in section III of this report) and Rahman et al. (1989c) 
(described in section V of this report). These studies considered both 
supercritical and subcritical flow regimes and used numerical techniques to 
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develop computational procedures for thin film flows. 

Thomas et al. (1989) approximated the flow to be one- dimensional m 
nature by neglecting any variation of velocity across the thickness of the 
filn. The continuity and momentum equations .ere integrated across the 
thin film to develop a single equation for the film velocity. The 
resistance to the flo. due to .all shear stress .as expressed in terms of a 
friction coefficient. Suitable expressions for the friction coefficient 
.ere developed for regions close to the entrance .here the thickness of the 
.omentum boundary layer is small compared to film thickness, and for 
regions far a.ay from the entrance .here viscous effects propagate all of 
the .ay across the film. These regions .ere matched .here the t.o friction 
coefficients yielded the same results. The equation of motion .as solved 
numerically using the MacCormack predictor- corrector method. Results .ere 
obtained for both stationary and rotating disks for a number of inlet 
Reynolds, Froude and Rossby numbers. In the presence of the jump, the 

outlet Fronde number .as al.ays assumed to be unity to simulate a situation 
.here the fluid experiences a free fall over the edge of the disk due to 
gravity. It .as found that jump moves do.nstrea. and may get .ashed a.ay 
.ith an increase in the flo. rate, angular velocity and inlet Froude 

number . 


The studies by Rahman et al. (1989a and 1989c) solved the flo. field 

numerically using a boundary- fitted coordinate system. No assumption 
regarding the velocity variation or friction coefficient .ere required. 
The studies, however, covered only the case of a stationary disk. Since 
the height of the free surface is dependent on the flo. conditions and is 
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not known ahead of tine, an iterative procedure had to be used to deternme 
the correct location of the free surface. The method presented by Rahman 
et al. (1989a) in section III of this report is termed as the ’pressure 

optimization method’. In this method, the shape of the free surface was 

assumed to be described by an algebraic equation with two or more arbitrary 
constants. The constants were optimized using an exhaustive search 
technique that minimized the difference between the computed free surface 
and ambient pressures. The method was found to be very satisfactory for 
zero-gravity flows, which was the primary emphasis in that paper. However, 
in the presence of gravity when a jump existed in the flow field, the 
method required the supercritical and subcritical flows to be computed 
separately, and the solutions were matched at the location of the jump. 
Even though this procedure yielded correct results in regions away from the 
jump, it could not account for the details of the flow field in the 
vicinity of the jump. This drawback was removed in the subsequent study by 

Rahman et al. (1989c). 

The method developed by Rahman et al. (1989c) in section V of this 
report is termed as the ’porous wall method’. In this method, the free 
surface was assumed to be a permeable wall through which fluid particles 
may leave or enter the control volume depending on the difference between 
the fluid and ambient pressures. The shape of the surface was corrected in 
successive iterations until the free surface conformed to a streamline and 
the penetration through the surface reduced to zero. This method computed 
the whole flow field as a single domain preserving the details of the flow 
in the jump region. The computed flow field, friction and heat transfer 
coefficients were reported in that paper. 


252 



The purpose of the present investigation Is twofold. First, to 
compute the flow field using the methods developed by Thomas et al. (1989) 
and Rahman et al. (1989a and 1989c) for the experimental conditions 
reported in Section II of this report. The theoretical and experimental 
results can then be compared to find the merits and drawbacks of the 
different prediction methods. The second objective of this work is to 
improve the existing theoretical models and to suggest a systematic 

computational procedure for thin film flows adjacent to a stationary or 

rotating disk. 

7.3 E flBiTIflNS OF MOTION 

The body-fitted coordinate system used for the two-dimensional 

numerical computation is shown in Fig. 7.1. The local coordinates are 

directed along lines connecting the centers of the adjacent grid cells. 

The z-axis is directed along the direction of radial spread, i.e., the main 
direction of flow. The y-axis is directed across the thickness of the 
film. The height of the free surface from the solid wall is denoted by 6, 
which varies with radial location on the plate. 

The equations governing the conservation of mass and momentum m a 

thin liquid layer involving a Newtonian, constant- property fluid can be 

written as 


dv = . i Vp + v r v + g 

ITt p y 
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Fig. 7.1. The coordinate s 
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The boundary 


conditions are given by: 


at entrance: 



(7.3) 


at exit: 


V - V 8=5 , i for stationary disk 

v ' out ’ out’ 

-+ 

- o, for rotating disk 
on 

b p = pg ( s ■ y) 


(7.4) 


at solid wall: 


, 0, for stationary disk 

l V , for rotating disk 
w’ 


(7.5) 


at free surface: ^ - 0> p - 0 

Here a is the coordinate directed normal to the surface under 
consideration. The p = 0 condition at the free surface is arrived at from 
a balance of normal stresses where the effects of surface tension and other 
stresses except for pressure is neglected. This is valid for a thin film 
flow with large values of the Reynolds and Veber numbers. 


7.4 Flow on a Stationary Disk 

The situation where the disk is stationary with fluid emanating at the 
center and spreading out umformly in the radial direction is considered 
first. Four sets of experimental data with flow rates ranging from < to lo 
lpm were taken where the film height was measured in the supercritical and 
subcritical regions. The measured data for 7 1pm and 11 1pm are plotted in 
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Figs. 7.2 and 7.3, respectively, vhere it is compared with the one- 
dimensional numerical solution by Thomas et al. (1989). It can be noticed 
that the experimental data compare reasonably well with the numerical 
solution in the supercritical region, but does not compare well in the 
subcritical region. The height of the jump is also not correctly predicted. 
In the numerical algorithm proposed by Thomas et al. (1989), the outflow 
fro. the disk was simulated by assigning a Tr = 1 condition at that 
location. Vhen the fluid flows out of the disk and experiences a free fall, 
the subcritical flow has to transform to supercritical flow, and a critical 

condition should exist somewhere in the vicinity of the outlet. However, 

assigning Fr = 1 at the outlet does not seem to be adequate in this case. 
From the figure, it can be also noticed that the measured height is higher 
than the critical height, which suggests that the film is held up by some 
sort of external force. It can be recalled here that in all previous 
computations, the surface tension was neglected in all regions of the flow, 

which may not be correct near the exit since an appreciable curvature of 

the free surface is encountered due to the rapid turning of the flow. 

In the presence of surface tension, the fluid and ambient pressures 
can be related by the equation 

_ 2<r (7.7) 

p e P a R 

where R is the radius of curvature. If the pressure is assumed to be 
hydrostatic in nature, this equation can be expressed in terms of film 

height, as 
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(7.8) 


where S is the actual film height at the exit, and S is the critical 
height predicted by the Fr = 1 condition at the exit. To solve this 
equation, an estimate for 8 is required, since that is unknown in general. 
In the case of the stationary disk, we may assume that the thickness of the 
film remains approximately constant in the region where the flow changes 
its direction from horizontal to vertical. This means that the radius of 

curvature remains of the same order as the film height. Substituting 

6 and solving for it, we obtain 
e 


* 

<5 + 


{ e = 


7 1 ST 

j 1 R 
1 


where 8 


*■[- 


2,r out’ /T 


2/3 


(7.9) 


equation (7.9), even though quite approximate in nature, seems to give 
a good estimate of film height at the exit. This is demonstrated in Table 
7.1, where the height calculated by equation (7.9) is compared with the 
measured height at the last data point from the exit where measurement was 
reliable. The deviation is within !«. From Table 7.1, it can also be 
noticed that the deviation has a sign change and increases with the flow 
rate. Vhen the flow rate is increased, the fluid particles move over the 
disk with a larger amount of inertia and get swept away further before 
experiencing the falling motion due to gravity. This increases the radius 
of curvature of the free surface and thereby decreases the effects of 
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Table 7.1: Comparison of Calculated and Measured Film 

Height at Outlet for a Stationary • Disk 


Flow rate Measured Height 

(lpm) (“O 


7 

9 

11 

13 


4.4196 

4.318 

4.369 

4.216 


Hal ml at ed Height 
(mm) 

4.3378 

4.4357 

4.5284 

4.617 


Deviation 

(70 

- 1.85 
2.73 
3.65 
9.51 


260 



surface tension. So, using R 
exit. The opposite is true for 


5 overestimates the film height near the 
e 

a smaller flow rate. The estimate, 


however, 


is reasonable considering the simplicity of the formulation 


The computed result with 6 - & e at the exit with the same inlet 

condition as the experiment is shown in Figs. 7.2 and 7.3 for flow rates of 
7 and 11 lpm, respectively, how, the subcritical film height compares 
reasonably well with the experimental data. In the supercritical regron, 
the height remains the same as that computed by the original procedure of 
Thomas et al. (1989). The location of the jump moves further upstream for 
both flow rates. The measured location of the jump is found to be somewhere 
in between the two numerical results. It can also be noticed that 
one-dimensional solution simulates the jump as a sudden rise of film height 
that occurs across a single grid cell. In reality, the jump occurs 
somewhat more gradually, which in our experiment spanned over a length of 
6-25 mm. Even with these drawbacks, the one- dimensional numerical solution 
developed by Thomas et al. (1989) may be used for a quick estimate of the 
subcritical and supercritical film height distribution if the modification 
of the outlet boundary condition proposed here is incorporated. 

The experimental data for 7, 9, and 11 lp» are compared with the 
two-dimensional numerical solutions of Rahman et al. (1989a and 1989c) m 
Figs. 7.4, 7.5, and 7.6, respectively. Here, we have used the first and 

last data point as the inlet and outlet conditions for the numerical 
computation. The computed results by the method of Thomas et al. (1989) 

are also shown here for a relat ive unison. In the two-dimensional 

numerical solution by Rahman et al. (1989a, 1989c), the velocity profile at 
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Radial distance, mm 

Fig. 7.4 Film height distribution over a stationary disk at 

flow rate of 7 Ipm 






Radial distance, mm 

Fig. 7.6 Film height distribution over a stationary disk at 

flow rate of 11 Ipm 




the entrance is assumed to be parabolic in nature with a maximum at the 
free surface. This condition is not present in our experiment since the 
flow comes through a slot before being discharged onto the disk. So, at 
the entrance the maximum is expected to be somewhere halfway between the 
solid walls. After a short distance from the entrance, the flow is 
expected to evolve to a parabolic profile with the maximum velocity at the 
free surface, since the no- slip condition exists on the solid wall and the 
zero- shear condition is present on the free surface. Therefore, using the 
first measured data point as the entrance instead of the physical inlet is 
more justified for comparison with the two-dimensional numerical solution. 
Also, analogous to the one- dimensional solution of Thomas et al. (1989), 
Rahman et al. (1989a and 1989c) used Fr = 1 to simulate the outflow 
condition from the disk. As we have already discussed, this is erroneous 
and the height (or Froude number) at the exit needs to be estimated taking 
into account the surface tension effects. Equation (7.9) can be applied 
for this purpose. Since the effects of the surface tension at the outflow 
location have been already demonstrated in Figs. 7.1- 7.2, in these 
computations the actual experimental measurement has also been used for the 
exit boundary to eliminate any error associated with the theoretical 

estimate. 

It can be seen that the ’pressure optimization method’ developed by 
Rahman et al. (1989a) predicts the correct location of the jump. The film 
height in the supercritical and subcritical regimes are also very well 
predicted. In this method, the supercritical and subcritical flows are 
computed separately by strictly imposing the inlet and outlet boundary 
conditions, respectively. The solutions are then matched at the jump 
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location by preserving the conservation of mass and momentum at the jump. 
In the original method, the jump happens instantly across one grid cell. 
In the solution presented here, we have modified it to allow the jump to 
form across a number of cells that result in a minimum overall pressure 
difference between the fluid and the ambient. The location for the 
initiation of the jump was kept the same and only cells at the downstream 
side were added to it. The quantity optimized was still the RMS error m 
the free surface pressure normalized by the initial total head (see Rahman 
et al., 1989a). In this process, both the supercritical and subcritical 
regimes were considered together along with the jump region. The solutions 
plotted in Figs. 7.4- 7.6 correspond to the slope of the jump that resulted 
in the minimum RMS error of the free surface pressure. The figures 
indicate that the two-dimensional solution predicts the location and slope 
of the jump somewhat better than the corresponding one- dimensional 
solution. So far as the location of the jump is concerned, the 
one- dimensional solution appears to be good for 7 1pm, but it predicts the 
jump at a radial location closer to the inlet for both 9 and 11 1pm. 

Even though the modified ’pressure optimization method’ gives the 
gross structure of the jump, the details of the free surface in the jump 
region can only be obtained by using the ’porous wall method’ developed by 
Rahman et al. (1989c). In this method, the local free surface height is 
changed in each iteration depending on the local velocities at the free 
surface. Moreover, a global variation is taken into account by minimizing 
the overall rate of penetration through the free surface. The surface 
evolves by itself as the iteration proceeds and the solution corresponding 
to the minimum rate of penetration is taken to be final. It can be noted 
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in Figs. 7.4- 7.6 that the film height predicted by the ’porous wall method’ 
most closely matches the experimental data, particularly in the jump 
region. As mentioned by Rahman et al. (1989c), the ’porous wall method is 
quite sensitive to the initial guess of the free surface height 
distribution and has rather slow convergence characteristics. So, the best 
prediction approach one can follow is to use the pressure optimization 
method’ first to get the gross structure of the jump, and refine that 
prediction using the ’porous wall method’. 

The details of the flow field for the flow rate of 11 1pm is shown in 
Fig. 7.7. It can be noticed that the velocity profile is parabolic at 
r = 0.108 m. This location corresponds to the supercritical region. Here 
the flow is well- structured with the maximum velocity at the free surface. 
This is the typical velocity profile in most of the supercritical region. 
However, as the jump is approached and the film height increases, the free 
surface tends to experience more and more force due to the gravitational 
head, and the location of the maximum velocity moves down from the free 
surface to somewhat inside the film. This can be seen in the vector plot 
centered around r = 0.153 m, where the jump has already started. After 
initiation of the jump the flow also loses inertia rapidly due to the 
sudden change in film height. Flow separation is encountered both at the 
solid wall and on the free surface, which can be noticed in both plots 
corresponding to r = 0.153 and r = 0.167 m. The separation from the solid 
wall is due to frictional resistance and was previously seen in the 
experimental work of Nakovyakov et al. (1978) and the numerical computation 
of Rahman et al . (1989c). Nakovyakov (1978) actually measured the wall 
shear stress for a circular hydraulic jump of a thin liquid film and found 
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3.0000E-01 m/s 


r=0.108m 



2.5000E-01 rn/s. 


r=0.153m 



2.2222E-01 m/s 


r= 0.167m 


Fig. 7.7 Velocity vectors for thin film flow adjacent to a stationary disk 
(flow rate — 11 1pm) 
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it to be negative for a short distance downstream from the location of the 
jump. Rahman et al. (1989c) reported a numerically computed friction 
coefficient that was negative in the vicinity of the jump. The flow 
reattaches to the surface after some distance in the subcntical region. 
The separation of the flow at the free surface is due to the gravitational 
pressure gradient and is more commonly termed as a ’surface roller in the 
hydraulics literature (Rajaratnam (1967) , for example) . This roller was 
not observed previously in other investigations on thin film flows. The 
reason was that a jump in a thin film flow was believed to be an 
instantaneous phenomenon and was modeled accordingly, so the details of the 
flow structure at the jump could not be captured. This phenomenon can be 
seen when the jump is allowed to develop gradually over a distance. In 
Fig. 7.7, it can also be noticed that the surface roller starts up at an 
earlier location and covers a larger region than the separation eddy near 
the solid wall. As seen in Fig. 7.6, the length of the subcritical region 
for this flow rate is quite small. If sufficient length is allowed, the 
flow again develops to a regular parabolic structure with the zero velocity 
at the solid wall and the maximum velocity at the free surface. 

7.5 Flow on a Rotating Disk 

When the disk rotates about its axis, in addition to inertial, 
gravitational and viscous forces, the centrifugal and Coriolis forces shape 
the flow structure as it moves downstream. For an appreciable rotational 
velocity (as is considered in our experiment), the subcritical flow region 
is swept out from the disk and no jump can be present. It may be recalled 
here that the length of the subcritical flow regime was small even with a 
stationary disk. If a disk of much larger diameter was considered, a jump 
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could be possible at a smaller rate of rotation. 


Ve first consider the data corresponding to an intermediate flow rate 
and angular velocity (11 1pm, 100 rpm) , which are plotted in Fig. 7.8 along 
with the theoretical prediction of Thomas et al. (1989). It can be noticed 
that the numerically determined film height is reasonable at smaller radii, 
but seems to deviate from the measurement at larger radii. Both the 
experimental measurements and the theoretical computation show that the 
film height increases downstream, attains a peak and then decreases. The 
experimental data also shows an increment of height near the exit. This 
increment is believed to be caused by the surface tension as discussed 
earlier in this section. The effect of surface tension is smaller for a 
rotating disk than that for a stationary disk since the flow moves over the 
disk with a higher velocity and results in a larger radius of curvature. 
In the numerical solution of Thomas et al. (1989), the flow was assumed to 
be strictly radial in nature with a superimposed solid- body rotation. The 
resistance to the flow due to friction was taken to be the resistance in 
the radial direction. In the azimuthal direction, the velocity remained 
constant all across the thickness of the film with no resistance from the 
solid wall. In reality, however, the velocity is expected to change due to 
the finite viscosity of the fluid, so there will be frictional resistance 
from the wall in the azimuthal direction. 

To account for the frictional resistance due to the angular velocity, 
we may define the total shear stress as 
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Radial distance, mm 

Fig. 7.8 Comparison of experiment with numerical solution of 
rotating disk (flow rate = 11 Ipm, rpm = 100) 




i 


(7.10) 


T = 
W 


T <t> 


+ r 


Here r r = c f (1/2 p V 2 ) , where c f can be calculated by using the Blasius or 
parabolic solution given by Thomas et al. (1989). can be estimated from 
the exact solution of laminar flow adjacent to a rotating disk in an 
infinite extent of fluid. As given by Schlichting (1979), 


r, = 0.6 (tr Pi 2 J ! 2 ( 7 - 11 ) 

Using r y as the shear stress at the- solid wall in the formulation by Thomas 
et al. (1989), the film height distribution was calculated and is shown in 
Fig. 7.8. The predicted film height compares reasonably well with the 
experimental data. So, at an intermediate flow rate, both components of 
frictional resistance are important and one cannot be entirely neglected in 
favor of the other. 


Figures 7.9-7.11 show the comparison of the numerical results with the 
experimental data for other rotational speeds tested in the experiment. 
For a smaller rate of rotation (55 rpm) , it can be seen that the modified 
numerical solution predicts the experimental data except for the region 
close to the outlet where surface tension is large. The original procedure 
of Thomas et al. (1989) somewhat underestimates the film height in the 
region away from the center where centrifugal force becomes large. For 
large rotational rates (200 rpm and 300 rpm), however, the method of Thomas 
et al. (1989) appears to be reasonably satisfactory and predicts the trend 
of the experimental data better. The assumption of solid- body rotation 
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Radial distance, mm 

Fig. 7.10 Comparison of experiment with numerical solution of 
rotating disk (flow rate = 13 Ipm, rpm = 200) 
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0 60 80 100 120 140 160 180 200 : 

Radial distance, mm 

Fig. 7.11 Comparison of experiment with numerical solution of 
rotating disk (flow rate = 7 Ipm, rpm = 300) 


inherent in the formulation of Thomas et al. {1989) become more appropriate 
at smaller values of the Rossby number, i.e., at higher rates of rotation, 
so the trend seen here may be expected. However, we should also keep in 
mind the limitations of a simple one- dimensional numerical solution. As 
seen in the experiment, waves appear on the surface and become more intense 
at a higher spin rates. The turbulence may also be triggered at a higher 
rate of flow or rotation. A two-dimensional solution with rotation may 
address some of these issues, which is beyond the scope of the present 
paper. 

7.6 CONCLUSIONS 

The experimental data were compared with the existing numerical 
solutions of radially spreading thin film flows adjacent to a stationary or 
rotating disk. Some of the drawbacks of the existing numerical methods are 
pointed out and corrected. The detailed structure of the flow for a 
stationary disk is also computed and discussed. 

It was found that surface tension is important at the outer edge of 
the disk. For a stationary disk, a simple model for estimating film height 
at this location is proposed. At other locations in the flow, the surface 
tension can be neglected. For a stationary disk, the one- dimensional 
numerical method of Thomas et al. (1989) or the two-dimensional methods of 
Rahman et al. (1989a and 1989c) should always be applied with the correctly 
estimated outlet Froude number, instead of using Fr = 1 to simulate a 
critical (free falling) outflow condition. 

In the case of a rotating disk, the frictional resistance due to 
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angular velocity was found to be important, it a smaller rate of rotation, 
the resultant frictional resistance was found to be more adequate to 
predict the film height than using only the radial component. 


The increment in film height at the jump was found to be gradual 

spanning over a number of grid cells. In the tno-dimensional -pressure 
optimization method’ after determining the location of the jump, an 
additional procedure for adjusting the slope of the jump may be included to 
arrive at a better structure of the free surface. The two-dimensional 


■porous wall method’ given by Rahman et al. (1989b) can be used as it is 
except for imposing the correct Fronde number at the outlet. It was also 
found that recirculating flow is present in the vicinity of the jump, both 

at the solid wall and on the free surface. The separation of flow at the 

solid wall is due to the viscous resistance exerted by the wall that 
overcomes inertia when the film height rises. The separation on the free 
surface is due to the adverse gravitational pressure gradient. 
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Appendix A: EQUIPMENT SPE CIFICATIONS 

* 

The equipment for the rotating disk unit is divided into 
subsystems for the description of the components: 

1. Rotating disk assembly 

2. Liquid film height measurement system 

3. Heat transfer measurement system 

4. Liquid circulation system 

1. Rotating disk assembly 

a) Rotating disk 

Material: 6061 aluminum 
Diameter: 406.4 mm 
Thickness: 6.35 mm 

- 4 

Surface finish: 3.81 x 10 mm 

- 3 

Perpendicularity to centerline of spindle: 2.5 x 10 mm 

- 3 

Flatness: 2.5 x 10 mm 

b) Collar 

Material: 316 stainless steel 
Outer diameter: 102 mm 

-4 

Perpendicularity to centerline of spindle: 6.0 x 10 mm 


four 
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c) Shim 

I 

Material: 316 stainless steel shim stock 
Thickness: 0.673 mm 

d) Spindle and motor 

• Vhitnon model 699-0080-000 "piggyback” spindle/motor assembly 
Spindle bearings: ABEC #7 angular contact ball bearings 

Spindle shaft: 440 stainless steel hardened to 50-55 Rockwell "C" 

except at shaft extension which will be 38-42 Rockwell "C" due to 
thin cross section 

Electric motor: 0.5 HP at 900 RPM with an input of 230/460 V, 3 

phase, 60 Hz (0.75 A @ 460 V) 

Driver and driven pulleys are interchanged to obtain speed range 
change. A frequency inverter will be used to vary the speed within 
each speed range. 

Frequency inverter: 4.5 A @ 460 V 
Low speed range: 55 to 374 RPM 
High speed range: 423 to 1987 RPM 

e) Radial encoder and speedometer 

• Airpax solid 60- tooth gear, magnetic pickup, and digital display 
Accuracy: ±1 RPM over speed range of 55 to 1987 RPM 
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2. Liquid film height measurement system 

a) Non- contact capacitance sensor 

• MTI model ASP-100-CTA 
Accuracy: ±0.002 mm 
Resolution: ±0.002 mm 

b) Data logging system 

• Fluke model 2280A Data Logger 

c) Digital micrometer 

• Microcode model 9598 
Accuracy: ±0.002 mm 
Resolution: ±0.002 ram 

d) Precision slide 

• Microslides model A-6166-LC crossed roller bearing slide assembly 

. 3 

Accuracy: Run-out over 130 mm = ±2.6 x 10 mm 

3. Beat transfer measurement system 

a) Etched foil heater 

• Minco mica heater 
Outer diameter: 356 mm 
Inner diameter: 102 mm 
Thickness: 0.635 mm 

(6) - 3.18 mm holes for thermocouples 
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Maximum operating temperature: 100°C „ 

Electrical specifications: 500 V, 41.67 fl, 12 A, 6000 V 

o 

Maximum heat flux: 6.6 V/cm 

Ceramic insulation 

• Aremco 502-600 machinable ceramic 
Maximum operating temperature: 590°C 
Thermal conductivity: 0.15 V/(m-K) 

Thermocouples 

• Type "T" copper- constantan 
Maximum useful range: -200 to 350°C 
Maximum error over range: 0.5°C 

Digital thermometer 

• Fluke model 2280A Data Logger 
Thermocouple conformity: Better than 0.05°C 

Slip- ring capsule 

• Electro- Tec model ETC P/N 67584 
Number of rings/circuits: 20 
Current ratings: 3 A at 500 VDC 

Contact resistance (noise): 25 mil at 100 mA (10 RPM) 
Life: 18 million revolutions at 8 RPM 
Operating environment: -55 to 125°C 


4. Liquid circulation systea 


a) Rotating union 

• Deublin series 1106 

b) Digital flow meter 

• Cole- Parmer model J- 5618- 04 licet semi- micro flow sensor 
Range: 1.51 to 17.0 LPM 

Maximum pressure: 3.76 x 10^ Pa at 105°C 
Linearity: *1 7 , full scale 
Repeatability: ±0.25 7 % full scale 

• Cole- Parmer model J- 5622- 35 Line powered flow meter 
Accuracy: ±1 7. full scale 

Repeatability: ±0.5 7. full scale 


286 



Appendix B: 


DESIGN OF THE ROTATING DISK ASSEMBLY 


The details of the design of the rotating disk assembly will be 
explained. This includes the rotating disk, spindle, electric motor, 
rotating fluid coupling, and the electric power supply for the motor and 

the etched foil heater. 


Figure B.l is an overall view of the rotating disk assembly as it appears 
in the experimental setup. The electric motor is mounted to the spindle, 
which in turn is mounted onto the vertical frame plate. The spindle is 
driven by a poly- groove belt which can be removed by taking a cotter key 
out and driving a pin down and out of a hole on the plate the electric 
motor is mounted on. This allows the electric motor to swing toward the 
spindle. This procedure is also performed to change the speed range of the 

spindle. 


Near the pulley on the spindle in Fig. B.l is the magnetic pickup for 
the tachometer. It is mounted through a housing such that it is 
approximately 0.254 mm (0.010 in.) from a special gear. As the teeth pass 
by the active sensor, a square wave proportional to the angular speed is 
generated, which is then read and displayed by the digital readout. 


The rotating disk and collar are also shown in Pig. B.l. The collar 
is mounted onto the spindle shaft with four screws which pass up through a 
ceramic insulator ring. This is to prevent excessive heat from reaching 
the high-precision bearings. An 0-ring seal is placed between the spindle 
shaft and the inner diameter of the collar to prevent water from passing 
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over the top of the collar. A shim is placed between the collar and 
rotating disk to set the initial gap height. Caution must be used when 
installing this shim so that the chamfer or bevel on the inner diameter of 
the shim is directed upward to mate with the small radius where the spindle 
shaft has a step in diameter. Otherwise, the shim and rotating disk will 
not seat properly. Also, it was found that the thickness of the shim does 
not correspond to the initial gap height, as was originally planned. 
Therefore, feeler gages must be used to directly measure this gap. Like 
the collar, the rotating disk has an 0-ring to prevent water from running 
directly down the spindle shaft. The large nut below the rotating disk 
holds the disk against the spindle. In order to achieve repeatability in 
the stationary film thickness measurements, marks were placed on the collar 
and the rotating disk for alignment purposes. The two marks were first 
aligned with each other, and then the marks were aligned with the 
capacitance probe. With this procedure, repeatability could be achieved. 

Some preliminary tests were made with two different collars which are 
shown in Fig. B.2. The first design was a reducer as shown in Fig. B.2a. 
The part of the reducer that was closest and parallel to the disk acted as 
a flow straightener, which directed most of the liquid to flow in the 
radial direction. However, part of the fluid would rebound off of the disk 
underneath the flow straightener, which resulted in splashing around the 
collar at moderate and high flow rates. To eliminate this problem, the 
collar shown in Fig. B.2b. was tested and it was found that the splashing 
was not present with the longer flow straightener. 


When the disk was rotated, it was found that a significant amount of 




Fig. B.2.a Collar with short flow straightener 



Fig. B.2.b Collar with long flow straightener 


Figure B.2. Collar design 
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water spray was generated by the liquid film leaving the disk and 
rebounding off of the Lexan water tank. This created problems with water 
droplets splashing onto the capacitance probe, so a screen mesh was draped 
on the sides of the tank to absorb the kinetic energy of the droplets and 
prevent the water spray from rebounding off of the tank. 

Below the rotating disk is a cylindrical coupling between the spindle 
shaft and the rotating fluid union. These items must be removed before the 
rotating disk can be taken off of the apparatus. 

Figure B.3 shows a view of the footprint where the spindle mounting 
holes are located. Also shown is the top view of the pulley housing on top 
of the rotating disk assembly. The low speed range is shown in Fig. B.3, 
whereas the pulleys are reversed to achieve the high speed range. 

Figure B.4 presents the side and top views of the aluminum rotating 
disk. Six holes are drilled from the underside to just below the surface 
of the disk for the placement of thermocouples along one radius. The top 
of the disk is coated with an aluminum oxide coating to prevent oxidation 
on the surface of the disk. 

Figure B.5 shows the design of the insulating ceramic disk. A cutout 
at the top of Fig. B.5 on the ceramic disk is for the electric heater 
leads. The six holes for thermocouples are drilled completely through the 
ceramic disk. 

The electrical diagram for the rotating disk unit is shown in Fig. 
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3.3. Side view #2 and top view of motor and spindle 
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I.’ lg ure B.6. Electric diagra m for the rotating distc unit 










B.6. The line voltage comes into the .receptacle to the left of the 
frequency inverter. The plug should never be disconnected from the 
receptacle when either of the transformers are operating. To energize the 
electric motor, first switch on the disconnect and then turn on the 
frequency inverter switch. Turn the system off by reversing the procedure. 
The electric motor speed is controlled with the remote control panel box 
supplied with the frequency inverter. After the start button is depressed, 
the speed dial can be turned to adjust the speed. The disk speed is 
leveled off when the "at speed" lamp is lit. 

To energize the electric heater on the rotating disk, turn down the 
variable ac transformer and turn on the main and heater disconnects. The 
variac can then be adjusted for the desired power input. Care must be 
taken because the output of the receptacle is 0 - 480 V. 
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Appendix C: ERROR AN ALYSIS . 


In order to properly report any experimental data, an analysis of the 
errors involved in taking the data must be given. In the present 
experimental setup, the errors involved in taking two measurements will be 
described and quantified: the liquid film height, and the liquid flow rate. 


Liquid film height 

The error in measuring the liquid film height involves three pieces of 
equipment: the digital micrometer, the capacitance sensor, and the 

rotating disk. To evaluate the error in the film height measurement, the 
individual errors of each of these pieces of equipment must be found. A 
summary of the method to determine the error will now be given. A full 
discussion on error analysis is given by Miller (1989). 

If N is a known function of n independent variables u^ , U 2 , ••• > u n > 

N = f(u 1 , u 2 , ... , u R ) 

the absolute error is given by expanding the function f in a Taylor senes 
f(u 1 ± Au 1? u 2 ± Au 2 , ... , % ± Au n ) = f(u r u 2 , ... , u n ) 


W[ Au l + ^2 AU 2 + ‘ 


+ d J- Au 
cm n 
n 


1 

5 


fit > 2 

75 ( lu i> 

a 


(C.l) 
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where the partial derivatives are evaluated, at the. measured values of u^ 
If the instrumentation errors Au^’s are small quantities, the terms which 
include (Au^ 2 will be negligible. Therefore, equation (C.l) is 
approximately given by 


f(u x * Au^ u 2 * Au 2 , ... , u n * Au n ) = f(u t , u 2 , 





df AlI 

V 2 


df . 


(C.2) 


The absolute error of the system is then given by 


E = AN = 
a 


Au l 




mi n 
n 


(C.3) 


where Au^ is the absolute error associated with the individual 
measurements. If the Au’s are considered to be statistical bounds on the 
instrument errors, such as 95 percent confidence levels, the formula for 
computing the overall error is modified using the root- sum square formula 


rss 


4 "l 


iU 2 


^ iU ” 


(C.4) 


The liquid film height is given by the following equation: 


6 = h’ + (h + Ah) - h" 


(C.5) 


where h’ is the digital micrometer reading, (h + Ah) is the capacitance 
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probe reading of the air gap between the probe and the dry disk, and h" is 
the capacitance p’robe reading of the air gap between the probe and the 
surface of the liquid film. 

Applying equation (C.4) gives 


E 


rss 


(Ah’)* 7 [A (h 7 Ah)] 7 T (Ah")^ 


(C.6) 


The first term in equation (C.6) is the stated accuracy of the digital 
micrometer, which was calibrated- at the factory with an NBS traceable 
Mark- Tech Laser Gage Model 7980. The maximum deviation over the range of 
movement reported on the calibration certificate is 0.001 mm. This value 
is then taken as the maximum error of the instrument. Therefore, 


Ah’ = 0.001 mm 

The second term involves the accuracy of the capacitance sensor, which 
is given by the manufacturer to be ±0.1 percent of the range when 
calibrated to a known standard. In the present experimental setup, the 
known standard was taken to be the digital micrometer. The capacitance 
probe was calibrated at the beginning of each test against the digital 
micrometer in the following manner: 

1. The capacitance probe is lowered until it lightly touches the surface 
of the disk, which can be seen with a lamp directed behind the probe. 
The capacitance probe is then zeroed at this point with the offset 
adj ustment . 
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2. The probe is then raised 0.0254 mm above the disk, where the digital 
micrometer is zeroed. This is done to eliminate any backlash in the 

digital micrometer. 

3. The capacitance sensor is then raised nntil it reads the maximum value 
of its range, which is 2.54 mm. 

The value read on the capacitance sensor readout is compared to what 
is shown on the digital micrometer readout. Any discrepancy can be 
eliminated by a gain adjustment which is supplied on the capacitance sensor 
readout. This sets the slope of the linear curve which relates the 
strength of the electric field to the height of the air gap. By choosing 
points at 0.0254 and 2.54 mm, the slope can be set over 99 percent of the 
range for the greatest accuracy. The linearity of the capacitance probe 
over its range was checked, and the maximum deviation was found to be 0.008 
mm. This value will be added to the stated accuracy of the probe. 

The second term in equation (C.6) takes on different values when the 
disk is stationary and when the disk is rotating. Vhen stationary, the 
second term is the accuracy of the capacitance sensor given by the 
manufacturer, which is *0.1 percent of the range when calibrated to a known 
standard, plus the maximum deviation of the linearity of the probe over its 
range. Therefore, when the disk is stationary, the second term is 

A(h + Ah) = 0.01 mm 

Vhen the disk is rotating, the error due to the disk not being 
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perfectly flat and the nobble caused by the bearings must be taken into 

account. The error when the dry disk is rotating is given by the sue of 

the standard deviation of the air gap measurement between the probe and the 
dry disk taken at the outer edge of the disk (where it is maximum) and the 
combined error of the capacitance sensor reported by the manufacturer and 
the deviation from linearity of the probe. The standard deviation of the 
disk with 50 samples is » = 0.02 Multiple populations with an 

increasing number of samples were taken to compare the means using a 
two-tailed a- test to ensure a sample number that is within the 95 percent 
confidence interval. Therefore, the maximum error in the tare data is 

A(h + Ah) = 0.03 nun 

The last term in equation (C.6) is again the error associated with the 
capacitance sensor: 

Ah" = 0.01 mm 

Due to the fact that the mean film thickness is being reported, 
further comments concerning the third term in equation (C.6) are necessary. 
At all flow rates and rotational speeds, waves were present on the free 
surface of the liquid film. Figures (C.l and C.2) present typical mean h" 
values and standard deviations as a function of the number of samples taken 
for the stationary subcritical and supercritical regions, respectively. In 
each case, the same population was examined for different numbers of 
samples to minimize the effect of the variation in the flow rate. For 
subcritical flow, a normal z-test was performed, and it was found that for 
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Figure C.l. Typical mean h M and standard deviation versus 
number of data points for subcritical flow 
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Number of data points 

Figure C.2. Typical mean h" and standard deviation versus 
b number of data points for supercritical flow 



125 and 250 samples, the means were equal at the 95 percent confidence 
level, so 250 samples were taken for each data point reported. Likewise, 
for supercritical flow 80 samples were taken in order to achieve the same 
confidence level. It should be noted that the standard deviations from the 
mean are not errors in the instrumentation but are actual variations in the 
mean film thickness due to the wavy nature of the free surface of the 
liquid film. A similar test performed when the disk was rotating showed 
that the number of samples needed for the 95 percent confidence level was 
the same as that for the stationary supercritical region. The overall 
root- sum- square error in the liquid, film height measurement is as follows: 

Stationary disk - 

Ah’ = 0.0001 mm, A(h + Ah) = 0.01 ram, Ah" = 0.01 mm 

E = ±0.01 mm 

rss 

Rotating disk - 

Ah’ = 0.0001 mm, A(h + Ah) = 0.03 mm, Ah" = 0.01 mm 

E = ±0.03 mm 

rss 

While the liquid film height at the collar could not be measured with 
the capacitance probe, the gap between the collar and the disk could be 
measured with feeler gages. Therefore, for the data point at the collar, 
the absolute error in the liquid film height was found to be 

E = ± 0 . 0 1 mm 
a 
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Flov «eter 


It is desired to determine the accuracy of the flow rate reported m 
the present experiment. The turbine flow sensor and digital flow meter 
were calibrated in the following manner: 

1. A line was drawn on the inside of a large bucket, and the volume to 
this line was measured with a 1000 * 5 ml graduated cylinder. The 
volume was measured three times and the average of these three 

readings was 17.23 ± 0.08 L. 

2. The flow rate through the flow sensor was allowed to reach the steady 

state at 12.5 ± 0.05 LPM, which is within the range of the reported 

flow rates. 

3. The timer was started and the water was directed into the calibrated 
bucket . 

4. The amount of time to fill the bucket to the line was noted and the 
actual flow rate was checked against the reading on the digital flow 

meter. 

5. The gain on the digital flow meter was changed until the flow rate 
read on the flow meter checked to within ±0.05 LPM of the measured 

flow rate. 

The accuracy of the actual flow rate measured with the bucket and the 
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stop watch can be determined with the method outlined in the previous 
section concerning the capacitance sensor. The equation for finding the 
flow rate is Q = V/t, where V is the volume and t is the amount of time. 
The absolute error is given by the following equation: 

E a = |!H + |H 4t | (c - 7) 

The time and volume measured are as follows: 

t = 1.37 ± 1/60 min. 

V = 17.23 ± 0.08 L 

The partial derivatives are 

£ - 07 ■ °- 730 lPI / 1 

B = - L = - ! 7 . 2 3 = _9 18 LPM/min 

91 t 2 (1.37) 2 

The error in the flow rate measured with the stop watch and the bucket is 
E = (0.730) (0.08) + (- 9.18) (1/60) = 0.2 LPM 

The maximum deviation of the flow rate read on the digital flow meter 
from the actual flow rate was found to be ± 0.05 LPM. The variation of the 
flow rate due to the circulation pump for each test where the film 
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thickness is measured was held to ± 0.05 LPM. Therefore, the total error 
on the flow rate is the sum of these three readings: 

E = ±0.3 LPM 
a 
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Appendix D: F.IPKRI1EWTAL DATA 


Experimental parameters: 


r . = 50.8 mm 

in 

r , = 203.2 mm 
out 

8 . = 0.267 mm 

10 o 

Fluid temperature = ambient temperature = 22 C 

Working fluid = de- ionized water 
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Table D.l . Experimental data for the stationary disk 


Radial Location 

7 LPM 

9 LPM 

11 LPI 

13 LPM 

15 LPM 

50.80 

0.267 

0.267 

0.267 

0.267 

0.267 

76.20 

0.322 

0.282 

0.269 

0.264 

0.250 

81.28 

0.330 

0.289 

0.277 

0.272 

0.254 

86.36 

0.345 

0.300 

0.284 

0.277 

0.259 

91.44 

0.361 

0.315 

0.297 

0.287 

0.266 

93.98 

0.378 

— 

— 

— 

" — ~ 

96.52 

— 

0.330 

0.310 

0.297 

0.267 

101.6 

— 

0.350 

0.325 

0.312 

0.288 

106.7 

— 

0.368 

0.340 

0.330 

0.301 

111.7 

— 

0.391 

0.361 

0.348 

0.320 

114.3 

— 

0.429 

— 

— 

— 

116.8 

— 

0.444 

0.383 

0.371 

0.340 

121.9 

— 

— 

0.411 

0.396 

0.359 

127.0 

4.37 

— 

0.442 

0.424 

0.388 

129.5 

4.42 

— 

— 

— 

... 

132.1 

4.42 

— 

0.490 

0.460 

0.412 

134.6 

4.47 

— 

— 

— 

... 

137.2 

4.47 

2.92 

— 

0.485 

0.441 

139.7 

— 

3.30 

— 

— 

— 

142.2 

4.47 

3.55 

— 

0.516 

0.469 

144.8 

— 

3.78 

— 

0.549 

— 
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Table D.l Experimental data for the stationary disk, continued 


Radial Location 

7 LPM 

9 LPM 

147.3 

4.52 

3.94 

152.4 

4.52 

4.14 

154.9 

— 

— 

157.5 

4.52 

4.24 

160.0 

— 

— 

162.5 

4.52 

4.34 

167.6 

4.52 

4.34 

170.2 

— 

— 

172.7 

4.52 

4.39 

175.3 

— 

— 

177.8 

4.52 

4.39 

182.8 

4.52 

4.39 

187.9 

4.47 

4.39 

190.5 

— 

4.34 

193.0 

4.42 

4.32 

195.5 

4.27 

— 


11 LPM 

13 LPM 

15 LPM 

— 

— 

0.487 

2.82 

— 

0.523 

3.20 

— 

— 

3.45 

— 

0.550 

3.63 

— 

— 

3.84 

— 

0.568 

4.09 

2.87 

0.596 

— 

3.17 

— 

4.19 

3.48 

0.644 

— 

3.63 

— 

4.29 

3.78 

— 

4.34 

4.04 

— 

4.37 

4.16 

— 

4.37 

4.22 

1.90 

4.32 

4.16 

2.22 

4.16 

4.06 

2.45 
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Table D.2. 


Experimental data for the rotating disk: 55 RPM 


Radial Location 

7 LPM 

9 LPM 

50.80 

0.267 

0.267 

76.20 

0.328 

0.317 

81.28 

0.333 

0.325 

86.36 

0.353 

0.335 

91.44 

0.366 

0.348 

96.52 

0.386 

0.363 

101.6 

0.404 

0.376 

106.7 

0.424 

0.396. 

111.7 

0.444 

0.409 

116.8 

0.462 

0.439 

121.9 

0.478 

0.452 

127.0 

0.490 

0.467 

132.1 

0.500 

0.477 

137.2 

0.513 

0.490 

142.2 

0.513 

0.513 

147.3 

0.513 

0.513 

152.4 

0.505 

0.518 

157.5 

0.505 

0.521 

162.5 

0.495 

0.533 

167.6 

0.480 

0.531 

172.7 

0.457 

0.508 

177.8 

0.457 

0.511 


11 LPM 

13 LPM 

15 LPM 

0.267 

0.267 

0.267 

0.279 

0.284 

0.274 

0.284 

0.287 

0.279 

0.292 

0.295 

0.287 

0.305 

0.302 

0.292 

0.312 

0.310 

0.300 

0.328 

0.325 

0.315 

0.343 

0.338 

0.330 

0.358 

0.351 

0.345 

0.373 

0.371 

0.363 

0.399 

0.396 

0.389 

0.414 

0.411 

0.401 

0.429 

0.444 

0.429 

0.455 

0.457 

0.442 

0.470 

0.483 

0.460 

0.500 

0.503 

0.475 

0.503 

0.518 

0.505 

0.528 

0.528 

0.518 

0.528 

0.543 

0.531 

0.528 

0.543 

0.541 

0.528 

0.546 

0.551 

0.528 

0.561 

0.551 
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Table D.2. Experimental data for the rotating disk: 55 RPM, continued 


Radial Location 

7 LPM 

9 LPM 

11 LPM 

13 LPM 

15 LPM 

182.8 

0.429 

0.508 

0.528 

0.546 

0.577 

187.9 

0.457 

0.508 

0.551 

0.546 

0.577 

190.5 

0.508 

0.551 

0.551 

0.561 

0.577 

193.0 

0.805 

0.729 

0.635 

0.587 

0.577 

195.5 

1.32 

1.09 

0.808 

— 

0.630 
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Table D.3. 


Experimental data for the rotating disk: 100 RPM 


Radial Location 

7 LPI 

9 LPM 

50.80 

0.267 

0.267 

76.20 

0.317 

0.312 

81.28 

0.325 

0.320 

86.36 

0.333 

0.328 

91.44 

0.348 

0.338 

96.52 

0.363 

0.351 

101.6 

0.381 

0.366 

106.7 

0.394 

0.378- 

111.7 

0.404 

0.394 

116.8 

0.414 

0.404 

121.9 

0.419 

0.417 

127.0 

0.424 

0.422 

132.1 

0.424 

0.427 

137.2 

0.419 

0.427 

142.2 

0.411 

0.427 

147.3 

0.401 

0.427 

152.4 

0.391 

0.424 

157.5 

0.383 

0.422 

162.5 

0.371 

0.411 

167.6 

0.358 

0.404 

172.7 

0.343 

0.399 

177.8 

0.335 

0.386 


11 LPI 

13 LPI 

15 LPI 

0.267 

0.267 

0.267 

0.300 

0.297 

0.307 

0.302 

0.302 

0.312 

0.307 

0.305 

0.320 

0.315 

0.312 

0.325 

0.325 

0.323 

0.340 

0.338 

0.333 

0.351 

0.351 

0.345 

0.366 

0.366 

0.361 

0.381 

0.381 

0.378 

0.404 

0.394 

0.394 

0.399 

0.406 

0.411 

0.419 

0.424 

0.427 

0.442 

0.424 

0.442 

0.460 

0.432 

0.442 

0.470 

0.432 

0.452 

0.477 

0.432 

0.460 

0.483 

0.429 

0.460 

0.493 

0.429 

0.452 

0.495 

0.419 

0.457 

0.498 

0.422 

0.452 

0.485 

0.414 

0.449 

0.483 
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Table D.3. Experimental data for the rotating disk: 100 RPM, continued 


Radial Location 

’ 7 LPM 

9 LPM 

11 LPM 

13 LPM 

15 LPM 

182.8 

0.328 

0.381 

0.396 

0.445 

0.483 

187.9 

0.315 

0.371 

0.394 

0.432 

0.462 

190.5 

0.312 

0.371 

0.409 

0.442 

0.460 

193.0 

0.310 

0.363 

0.406 

0.442 

0.465 

195.5 

0.330 

0.373 

0.406 

0.432 

0.457 
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Table D.4. Experimental data 


Radial Location 

7 LPM 

9 LPM 

50.80 

0.267 

0.267 

76.20 

0.302 

0.290 

81.28 

0.305 

0.297 

86.36 

0.307 

0.300 

91.44 

0.310 

0.305 

96.52 

0.307 

0.307 

101.6 

0.307 

0.312 

106.7 

0.305 

0.315 

111.7 

0.297 

0.315 

116.8 

0.290 

0.312 

121.9 

0.282 

0.305 

127.0 

0.272 

0.302 

132.1 

0.259 

0.292 

137.2 

0.249 

0.282 

142.2 

0.239 

0.272 

147.3 

0.229 

0.269 

152.4 

0.221 

0.254 

157.5 

0.211 

0.246 

162.5 

0.201 

0.231 

167.6 

0.193 

0.226 

172.7 

0.188 

0.216 

177.8 

0. 180 

0.213 


for the rotating disk: 200 RPM 


11 LPM 

13 LPM 

15 LPM 

0.267 

0.267 

0.267 

0.279 

0.269 

0.272 

0.284 

0.274 

0.279 

0.287 

0.277 

0.284 

0.295 

0.284 

0.292 

0.300 

0.287 

0.300 

0.305 

0.292 

0.310 

0.310 

0.302 

0.322 

0.315 

0.307 

0.330 

0.318 

0.318 

0.343 

0.318 

0.320 

0.350 

0.318 

0.328 

0.356 

0.315 

0.325 

0.358 

0.305 

0.320 

0.358 

0.302 

0.315 

0.356 

0.290 

0.307 

0.353 

0.281 

0.297 

0.348 

0.269 

0.295 

0.340 

0.259 

0.284 

0.333 

0.251 

0.277 

0.325 

0.241 

0.269 

0.320 

0.236 

0.259 

0.307 
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Table D.4. Experimental data for the rotating disk: 200 RPM, continued 


Radial Location 

7 LPM 

9 LPi 

182.8 

0.178 

0.203 

187.9 

0.173 

0.196 

190.5 

0.170 

0.196 

193.0 

0.180 

0.193 


11 LPM 

13 LPM 

15 LPM 

0.223 

0.244 

0.297 

0.223 

0.241 

0.287 

0.218 

0.236 

0.279 

0.213 

0.234 

0.279 
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Table D.5. Experimental data for the rotating disk: 300 RPM 


Radial Location 

7 LPM 

9 LPM 

11 LPI 

13 LPM 

15 LPM 

50.80 

0.267 

0.267 

0.267 

0.267 

0.267 

76.20 

0.305 

0.295 

0.289 

0.282 

0.281 

81.28 

0.302 

0.300 

0.289 

0.287 

0.282 

86.36 

0.295 

0.295 

0.287 

0.292 

0.285 

91.44 

0.292 

0.297 

0.289 

0.295 

0.289 

96.52 

0.289 

0.297 

0.292 

0.301 

0.295 

101.6 

0.284 

0.295 

0.289 

0.303 

0.302 

106.7 

0.279 

0.292- 

0.292 

0.305 

0.310 

111.7 

0.272 

0.287 

0.292 

0.304 

0.314 

116.8 

0.262 

0.282 

0.289 

0.305 

0.317 

121.9 

0.251 

0.274 

0.287 

0.303 

0.321 

127.0 

0.244 

0.264 

0.284 

0.303 

0.319 

132.1 

0.231 

0.259 

0.274 

0.297 

0.313 

137.2 

0.223 

0.246 

0.262 

0.288 

0.308 

142.2 

0.213 

0.236 

0.251 

0.280 

0.298 

147.3 

0.206 

0.229 

0.244 

0.279 

0.290 

152.4 

0.198 

0.218 

0.234 

0.264 

0.284 

157.5 

0.193 

0.203 

0.226 

0.250 

0.272 

162.5 

0.185 

0.201 

0.221 

0.250 

0.265 

167.6 

0.183 

0.195 

0.208 

0.239 

0.257 

172.7 

0.175 

0.185 

0.203 

0.230 

0.247 

177.8 

0.173 

0.178 

0.195 

0.227 

0.245 
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Table D.5. Experimental data for the rotating disk: 300 RPM, continued 


Radial Location 

7 LPX 

9 LPI 

11 LPM 

13 LPM 

15 LPM 

182.8 

0.165 

0.175 

0.193 

0.218 

0.233 

187.9 

0.160 

0.168 

0.183 

0.209 

0.219 

190.5 

0.160 

0.168 

0.183 

0.209 

0.227 


318 



